I am posting an updated crosslinking script here as the old thread is closed. Three updates are included as documented in the script header:
# JD 03/09/2015 Fix bug where maxbond turned off atom typing for ewald and atom based charges # JD 03/10/2015 Configured settings to run with or without GUI # JD 03/13/2015 On-the-fly conversion limit, eliminating problem of open epoxy rings from deleted xlinks
Here is the latest update dated "Sep 4 2015" which includes the following changes:
# JD 09/02/2015 System does not need to have molecule objects defined # JD 09/03/2015 Rename starting structure to avoid conflict # JD 09/03/2015 Count opt & md steps, de-fork analysis dynamics, new GUI option md steps per radius # JD 09/04/2015 Fix fragment analysis
I am pleased to release a new update to the crosslinking script. This one contains some significant new features as well as bug fixes.
The main updates are as follows:
The script no longer depends on molecule objects (see checkbox). This enables the simulation of certain interfacial systems as well as self-crosslinking silanes.
Support for polyurethanes
The changes are detailed in the script header: # JD 10/14/2015 GUI parameters are optional (allows simpler GUIs to be defined) # JD 02/01/2016 Allow run to continue even if molecules cannot be assigned (due to infinite network) # JD 02/03/2016 Enable single or multi condensation reactions for R1/R2 atoms with multiple OH groups # JD 02/10/2016 Prevent intra-molecular xlinks (multiple between same fragments) (incl concurrent) # JD 02/10/2016 Flag to run without using molecule objects (no stats output) eg if self-xlinking or interface # JD 02/29/2016 Switch intra-prevention method to pathlength basis rather than relying on molecule objects # JD 06/01/2016 Conversion formula with total R1 atom count on the fly in case R1 atoms are deleted during reactions # JD 06/02/2016 Global replacement of the term 'oligomer' for the more common term 'monomer' # JD 06/02/2016 Non-periodic systems (conditional pressure analysis), might work with 2D periodicity also # JD 07/12/2016 Fixed bug in self-crosslinking systems causing xlinking in excess of proper coord nr # JD 09/14/2016 Support for polyurethanes - convert C=N bond to single bond in isocyanate group # JD 09/15/2016 Fix bug with mol-based intra exclusions - appending index to xlinker mol names
Thank you very much for providing this helpful and nice crosslink script.
I am simulating the crosslinking procedure using the update script provided here.
I am used DGEBA as the monomer and TETA as the xlinker. The initial amorphous structure of monomer/xlinker was constructed by Accelrys amorphous builder in Materials Studio 8.0. Then I ran the script for the crosslink after changing the file names and making "FALSE" to "my $polyurethane". The script ran well until an error message occurred in line 1939 as following: Automatic charge group calculation failed. Atom does not belong to bond in Forcite.Energy (function/property "Run") at -e line 1939.
I have no idea about this error. Could you please give some comments on it?
Without some files to look at I cannot really say what is the issue. I can describe what the script is trying to do on line 1939. After new crosslinks are formed, the chemistry is adjusted and so the charge groups (if you are using charge groups) need to be reassigned. This is done via a Forcite energy run. You are getting this message because the charge group assignment failed.
You can take the structure and try a Forcite Energy run on it with the same settings. It should also fail. There might be some bad chemistry in there. Take a closer look at the crosslinks that have formed to see if there is anything amiss.
Since you are doing a standard epoxy-amine reaction, your chemistry settings should be like this:
$openRing = TRUE
$remove_condensation = FALSE
$polyurethane = FALSE
There may be some chemical cases in which charge groups cannot be calculated, in which case you could always switch to Atom based (less accurate) or Ewald (slower) charges. But I doubt it in this case as I have run many of these types of systems.
Thank you very much for your reply and detailed explanation of script in line 1939. Yes, it also failed in the Forcite Energy. I tried several times, checked the structure, and found a strange link as following:
The initial periodic structure was formed in Amorphous Cell module. After the crosslink in the first cutoff distance, the two atom within the period boundary were defined to be linked (R1-1 and R2-1 in red circles). But in the real structure they were not connected. I think this is the problem for the charge group calculation. If I try "Atom based", it sometimes works. However, I think the crosslinking structure is not right. Could you provide some suggestions on this problem?
As to another question I posted, I have chosen another forcefield and the structure of DGEBA was reasonable then. Thank you very much for your concern.
What forcefield are you using? That will definitely affect the charge group calculation. For DGEBA/TETA crosslinking I always use COMPASSII and it works very well. If possible, post the XSD file for the crosslinked structure above. It looks to me like there is a bond connecting R1-1 and R2-1.
I just try your Crosslink script in Materials Studio 7.0 using the sample input system DGEBFx64_TETAx32. I observe that there are many times of repeating the run in order to complete crosslinking. Can you give me an average real time which is characteristic for completing this procedure? Now, it is running about one hour and more and it hasn't finished yet.
I thought I responded a few days ago, but it looks like my post never made it. Anyways, I have looked at your structure and everything seems fine. The automatic charge group assignment using Divide and Conquer is successful. I am using Materials Studio 2017 and I think that is the difference here. Back in Materials Studio 8.0, I believe there was a problem was hence fixed and we saw it crop up in certain crosslink structures like these.
My advice then is to upgrade if possible. Otherwise, you will need to choose Atom based (faster) or Ewald (more accurate) for the charges.
Yes, there is a problem in assigning charge with Divide and Conquer. The process was good when running the initial radius (3), but the error occured when running the next radius (3.5, which I posted here). Then I tested the structure by auto-calculation in Charge Groups function with Divide and Conquer, and the error occured with informaiton "Atom does not belong to bond" in Materials Studio 8.0. I will try Atom based. Thank you very much for your comments and suggestions!
I would like to ask you if I can create structures of crosslinked chitosan and crosslinked hydrogel HA-TA (tyramine derivative of hyaluronic acid) using your crosslink script. Do you believe that I can't use it correctly due to my old version of Materials Studio 7.0? Also, I have an error message. I created a structure with Amorphous Cell using DGEBF and TETA monomers and when I call your crosslink script (last version Sept_2016), the error message is "Illegal division by zero at -e line 370" where is a calculation of conversion = 100*(reactedMonomerAtoms/totalMonomerAtoms). If you have any idea what it means I will be grateful.
The error message indicates that the script failed to find any monomer atoms in the system. It is looking for atoms whose name begins with R1. You can count them in the visualizer by using Edit | Atom Selection with Name matches R1*. Then look at the atom count shown in the Display Style dialog.
Regarding chitosan and HA-TA, I would have to look at the chemistry of the crosslinks to see if the script is applicable. Do you have an image of the specific crosslink structures of interest?
Moreover it is important to note that this type of script was really developed for thermoset polymers rather than elastomers or hydrogels. I would assume the latter have much lower crosslink densities which would require increasing the system size to possibly impractical levels. So I would recommend estimating the crosslink density and how many crosslinks you would end up with in a cell that is a few nm on a side.
Hi, Jason The script you provided is so helpful and it is really very selfless of you to do so. I really appreciate it. The simulation I gonna do is about crosslinks between Si-OH. Two "Si-OH" react and lose a "H2O", then "Si-O-Si". So I need to use the script for reference and write my own script. May I ask for some suggestions? Here are some questions I have considered. How to label the atoms when I create the input structure? Because every monomer has O and H which are reactive. And the close contact i need to consider is hydrogen bond, so i should take different way to make the crosslink operation.
Hi Jason thank you so much for your answer. it's very nice of you to respond to many other researchers as well as me. here i describe my situations more specifically. my simulation is about hydrogels and the monomer is something like "R-Si-(OH)3" or "Si-(OH)4". and the reaction has fundamental difference with that in your script, cause i have two reactive atoms in one monomer, not R1 in monomer and R2 in xlinker. and i also face the problem mentioned above by loannis, in line 370, illegal division. your answer reminds me of a mistake that i label O atom R1 and H atom R2. i realize that i was wrong. still, may i ask for further suggestions?
The script does support your case actually. In the script I refer to that as self-crosslinking since the R1 and R2 atoms are on the same molecule.
The correct way to run that system is to apply the following flags:
OpenRings = false
Condensation = true
Polyurethane = false
React_Mult_Monomer = true (where Si = R1)
React_Multi_Xlinker = false
Use_Molecules = false
For self-crosslinking systems it is critical to set Use_Molecules to false so it does not assume there are two separate species present. This skips over the analysis routines that are useful for tracking the statistics of ordinary two-component systems like epoxies.
The multi flag allows the Si (R1) atom to react multiple times (as many times as it has OH groups). There is one extra trick being used here. To prevent multiple crosslinks between the same two groups of atoms (intra), I am using a path length criterion to exclude such cases (which would otherwise be numerous in this type of system).
Let me know, and I will post some example Si(OH)3 systems that I have crosslinked.
Hi Jason Thank you so much for your detailed reply. That is what I do to change the script(Sep-15) according to your reply(maybe I didn't do it correctly): - my $noMol = FALSE - my $openRing = FALSE - my $remove-condensaion = TRUE - my $react-multi-monomer = TRUE - my $react-multi-xlinker = FALSE - my $polyurethane = FALSE And the bug is in the image.
And one more thing, I'm not a native English speaker. So when you said"Let me know", I'm confused. Let you know what?............or it is just a modal particle? sorry, i just can't get it.............
Hi Jason Please excuse me for keeping asking questions and your advice is really helpful. There are some new questions. If I want to start DPD simulation, how to define a bead when the polymer has crosslink structure? I mean crosslink structure is not like copolymer or homopolymer, so the defining way can be different. I guess if a bead consists of four(maybe 5, 6 or more)-ring (4 monomer)can be used. And my system is polymer[Si(OH)4 condensation] dissolved in water. So i think maybe i should define two types of beads which stand for polymer and water respectively? About the Flory-Huggins parameter, I know that we can run Forcite Module simulation to get "cohesive energy density"(let me say it Ec).
Is the solubility parameter equal to "sprt(Ec)"? I have learned the equations about how to calculate Chi parameter.
Thank you so much for your time.If my questions show my ignorance, really sorry for that.
Good question. But this thread is really for questions about using the crosslinking script itself. Could you please post you question under a new thread titled "Coarse-graining crosslinked structures" or something like that. The equation and your interpretation look fine, but lets take a closer look at how you plan to map atoms to beads given an arbitrary network structure.
Hi Jason I want do some modification to your crosslink script, specifically, the “sub XlinkSet” part. Following your advice, I successfully crosslinked "Si-OH" system, and the Si and O atoms are set as "crosslink". What I want to do is setting only "O" as "crosslink", and I failed rewriting the script several times. How can I do it?
Hi Jason This is Fang again After such a long time, I noticed that there was a mistake I might made. For a self-crosslinking system like "Si(OH)4", you said that it was critical to set "Use_Molecules" to False. And I can't find where the setting sentence in the script. I've successfully got the crosslink result, and it seems right. But I still want to make sure that I set it correctly. Thanks a lot Fang
well, I have a question though, in the "sampleInput" folder you create a file named "DGEBFx64_TETAx32", where in that you inserted 64 DGEBA and 32 TETA. that means there are 128 active sites on the DGEBA while 192 active sites on TETA. you didn't observe stoichiometric ratio in that mixture! Is there any reason for that?
By the way, I faced similar error by Anthea, I tried both atom based and group based for the charge and also applied COMPASS II forcefield. in the meanwhile in the real structure, the active sites were not connected, unlike the report by the script (just similar to Anthea case). could you attach a crosslinked structure of DGEBA/TETA or DGEBA/DEDTA system using this script? It would be very useful to check the structure and option you apply. by the way, I use Biovia Material studio 2017 (22.214.171.124)
You are right, 192 is the number of possible TETA crosslinks (6 each) if you check the React_Multi_Xlinker box. Otherwise each amine can only react once, which gives your 4 per TETA and a total of 128.
I have attached a 90% crosslinked structure of DGEBF and TETA. Could you describe the problem you are seeing a bit more? Or possibly attach an example?
Here is the latest update to the script. I have also updated the Usage Guide significantly.
# JD 12/29/2016 Fix bug causing over-crosslinking in self silanol Si(OH)3 systems # JD 12/30/2016 Prevent user from selecting Ewald for nonperiodic systems # JD 01/04/2017 On-the-fly conversion calculation that works even when R1 atoms are variable # JD 01/04/2017 User can specify ensemble (normally NPT)
Thank you for the attached file. I opened the structure and I found out similar ambiguity. please look at the structure in the circles. The active site on the TETA reacted to DGEBA, but in the real structure, it seems to have no connection!! am I right?
I enclose the final XSD file and the perl code. In the perl code, you can find the setting I used for the crosslinking procedure.
I use DGEBA as the monomer and DEDTA as the crosslinker (I attach both structure, too).
I also attach my starting cell which was built by amorphous cell construction and several NVT and NPT simulations.
The mixing ratio of DGEBA/DEDTA was 8:4.
Is there any problem with the final structure? it looks like that there are some bondings with no connection (Same as my previous problem in your structure)? using "select fragment" for the structure with fake bonding also prove my claim. when you select the DGEBA (for example) which shows fake crosslinking with xlinker, no fragment is selected! it shows that no crosslinking occures in that DGEBA, unlike the graphical observation.
I think the problem comes from the stepwise linking procedure in which per higher linking distances the structure couldn't reach its equilibrium bond length of the bonding due to movement restrications of the whole molecules and structures.
By the way, if we turn on the annealing step, it will perform after each crosslinking radious. is that right?
It looks fine. Both of the amines you circled are drawn with a half-bond. The other half of each bond can be found on the other side of the cell: Periodic boundary conditions.
As an aside, I have used a labeling scheme for the names of all crosslinked atoms. You can expect to find a crosslink between R1-n and R2-n, where n is substituted by a number indicating the crosslink index. That way it is easy to use the atom selection tool to search for crosslink pairs. For instanced searching for R?-1 will select the pair of atoms involved in the first crosslink that was formed.
The Select Fragment tool only applies to what is shown in the Display Range. If you open the Display Style dialog, you'll notice that the Lattice style is Original. This displays the atoms in so-called "diffusion coordinates". Thus it is possible to visually track how far they have moved relative to their original positions. Notice that the particular fragment you have focused on has moved the farthest from the cell, because it did not form any crosslinks until late in the game: #12 & #13.
Here are a few experiments you can try:
Increase the lattice display range to show multiple cells. i.e. decrease the min and increase the max values. Now try Select Fragment again. Notice that as the display range expands, your fragment selection gets ever larger. That is because, as indicated in the Progress.txt file, you have a single network fragment in this system, so it is infinite.
Go back to the starting view. Select the xlink bond named "xlink-12". It is one of the bonds attached to R2-6-12. Now we ask, where is the other atom that this bond is connect to? Let's find it. On the lattice display range, change A min to -1 and C max to 2. Without deselecting the bond, recenter the view and zoom in. Voila, the missing atom appears!
Go back to the starting view. Change the lattice display style to Default. Now everything is in the box. Use Edit | Atom Selection to select the atom named "R2-6-12". Recenter and zoom on this selection and notice that you can now see the full crosslinks, with the R1-12 and R1-6 atoms on the other sides.
Create an empty structure. Create a new box of the default size. Add a pair of atoms at fractional coords 0.1,0.5,0.5 and 0.9,0.5,0,5. Move them a little closer and calculate bonds. You will see the same effect with the half bonds on each side.
Thank you very much for the complete explanations. It was very very very .... very useful and eliminated any ambiguity.
As the last question, if we turn on Temperature cycle ($UseTempCycle = TRUE;), then this step will perform after each crosslinking radius? right? Do you see any difference in the final crosslinked structure with or without temperture cycle? Is it really useful to turn on temperature cycle?
I am glad it was helpful. For your question concerning the temperature cycle, that is something we tried out but I have never really found that it is worth the extra cpu time. I cannot find much difference in the networks being formed. I prefer to run at a high temperature 500-600K the entire time.
By the way, I noticed that last crosslink in your Progress.txt resulted in the maxBondEnergy of about 20 kcal/mol. That is a bit a high, indicating some stress. The maxbond energy should be below 10 normally. Try running some more dynamics to see if you can relieve the stress.
Dear Jason About the # JD 12/29/2016 Fix bug causing over-crosslinking in self silanol Si(OH)3 systems
I guess you were talking about over cross linking about Si(OH)4 or R-Si(OH)3? Indeed, in the final xsd of an old version crosslinking script job, there are Si atoms which are connected with 5 bonds. Thank you Jason, you really did great job!
Yes exactly. I have fixed the bug and also added a coordination check to identify such problems pro-actively in the future. Thus if you have atoms with a coordination number of 5 you will see a warning in the log file.
I think there is a bug of conversion calculation algorithm, when the system is about Si(OH)4. And it's just my personal opinion, maybe it's wrong.
Let's assume this situation: There are 100 Si(OH)4 in box, i.e. 100 Si and 400 O. And the crosslinking result shows that conversion is 80%. So there should be 80 Si atoms that are reacted. Due to the multiple reaction of Si, there will be more than 80 xlink bonds created. Actually, the number of xlink bonds are likely to be close to 200. Because every xlink bond consumes 2 -OH, the max number of xlink bond is 200. Let me assume the xlink bond number is 180, so there will be 180 O in the system which are reacted. And the number of Si atoms which are connnected to the 180 O atoms is definitely to be more than 80. Let me assume the number is 95. So the real result is, there are only 5 Si atoms which are absolutely unreacted. And there are 15 Si atoms which, although, is not involved in xlink reaction, and connected to a reacted O atom, i.e. there 15 atoms are in crosslinking strucure.
My point is, even though a Si atom is unreacted, it also may be connected to a O atom which is reacted. So the conversion result are very likely to be smaller that the real conversion. And we can't find those unreacted and connected Si atoms, although we can know the amout.
Thank you again for your great script and updates. I have two more questions:
1- could you please explain "use_molecules" option more clearly?
2- Are the two files "Presentation: Simulating Crosslinking in Materials Studio ...." and " New crosslinking feature in Materials Studio Collection" available in the net? if no, please kindly attach the files.
Hi Jason If R1 has two close contacts with R2-1 and R2-2 at the same time, and the distance are 3.1 Ångstrom and 3.2Ångstrom respectively, which will be reacted first? Thank you for your answering Fang
Use_Molecules: Check this box to use molecule objects to define the crosslinking entities in the system, thereby enabling a statistical analysis of the crosslinking through the course of the cure. Do not check this box if the R1 or R2 atoms are attached to infinite networks such as a substrate that would not allow molecule fragments to be defined. Also, do not check this box for self-crosslinking systems (i.e. where the monomer and xlinker are the same), such as R-Si(OH)3. In these cases, the script uses a different algorithm to identify suitable crosslink candidates based on the path length between R1 and R2.
In the original version of the script, it was assumed that R1 and R2 belong to distinct monomer and hardener molecules. Here you need to understand what Materials Studio considers to be a ‘molecule’. Molecule objects may or may not be defined within a given MS structure. To be defined, there needs to be a group of atoms that is linked together. However a molecule is not necessarily the same thing as a fragment. In fact, as the monomer and hardener molecules bind together and begin to form a network, the original molecule entities are unchanged.
This fact was convenient for a couple of purposes in the original crosslinking script. (1) We wanted to prevent ‘intramolecular’ crosslinks. This is the scenario where a pair of molecules has more than one crosslink between them. (2) We wanted to generate some statistics to indicate for instance, the mean number of crosslinks on each monomer. Thus the original script was built around the molecule concept.
Then of course some cases came along that violated our molecule assumption. (1) What if you are crosslinking some monomers to a surface? (2) What if the monomers have both R1 and R2 atoms and there are no hardener molecules. For these reasons, I creates a no-molecule option. To solve the intramolecular problem, I used a path length criterion instead. For the statistics, I just turned them off. So you can always use the no molecule option, you will just get less information. In the cases of surface or self crosslinking you must use it.
Let me address your question on conversion of May 13. The way I calculate conversion is the percentage of R1 atoms that has reacted. You are right that in some cases there could be more crosslinks than reacted R1 atoms, since R1 atoms can react multiple times (with the React_Multi_Monomer flag).
Let’s take a look at your example of 100 Si(OH)4 self-crosslinking molecules. I agree that the maximum number of crosslinks that can form is 200, since for each crosslink, one OH is lost via the condensation reaction. So if the simulation produces 180 crosslinks and a conversion of 80%, that simply means that you have an average of 180/80 = 2.25 crosslinks per Si atom. But that is impossible, since the average must be less than or equal to 2. Thus the maximum number of crosslinks that can form at 80% conversion is 160.
The above assumes of course Si=R1 and O=R2. It would not be such a good idea to do it the other way around since the conversion calculation would be based on a species with a declining population. I am not sure that there is a universal definition of crosslink conversion, but what I am aiming for here is something that is unambiguously defined and is a reasonable indicator of progress of the network formation. You are welcome to define your own metrics.