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.
I need to output the trajectory file, so I try to add something to your script. As far as I'm concerned, I need to output a ".xtd" everytime before the order of "$result -> Trajectory -> Delete". Is that right?