I can't figure out how to use restraints correctly, and would really appreciate some help. I simply want to keep some small molecules from getting close to the box boundaries.
What I've tried: I have a simple 31 x 31 x 31 A box filled with water and five hydronium molecules. I do mini abnr nstep 300 first and then successfully run a dyna CPT simulation without errors. Now I want to keep those five hydroniums from reaching the box boundary. So I tried to harmonically restrain them to the box center:
CONS HMCM FORCE 100.0 WEIG REFX 15.5 REFY 15.5 REFZ 15.5 - SELECT SEGID HYDR END
I placed this command in various positions in the input file: - anywhere in the dynamics block: I just get "extraneous characters found" in the output - after mini, before dyna: all molecules freeze except hydronium
I don't want the restraint included in the minimization, I only want it active during a regular MD run. I also played around with CONS HARM ABSO. I read the cons.doc but it doesn't tell me where to place the command.
Questions: 1) Where and how do I write the restraints block in the input file for a dynamics run? 2) Is this even a good way to keep my molecules in the box center? 3) Does RDIST of HMCM keep molecules within or outside a specific radius? 4) Is there a way to specify the box center instead of a reference coordinate set to restrain to with CONS HARM ABSO?
This seem to be basic question and I apologize if I overlooked something in the docs or the forum posts. I spent several days searching, but I am still a CHARMM beginner. Any advice or sample input will be greatly appreciated.
Hi, CHARMM executes each command as it is encountered in the input file, so your concept of "blocks" is somewhat off. The CONS command sets up the data structures required to apply the restraint, and all energy calculations aftermthe CONS command will include the restraint energy ( ENERgy, MINImize, DYNAmices) That said, if you want the restraints to be applied in the dynamics, but not in the minimization, you simply put the CONS command in between: MINImize .... CONS ,.. DYNAmics....
CHARMM allows you to apply several different restraints. Without knowing what your goal is, it is difficult to give adequate advice. Why do you want to keep your hydronium ions in the center? They cannot escape the box - you do apply PBC. And why do you play with hydrogen bond definitions (DONOR REMOVE; HBONDs...? There is NO hydrogen bond energy term n current CHARMM force fields.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
Many problems are resolved by carefully reviewing the output log file; you will need to learn how to do this if you wish to advance.
Besides what Lennart has already noted, I should point out a couple errors in the input:
The TCONS keyword conflicts with and overrides the HOOVER keyword, and invokes the older, flawed Berendsen thermostat; you should remove the TCONS keyword.
The system may have a net non-zero charge with Ewald, which is not recommended, esp. for a beginner, as QCOR 0.0 is probably not the best choice for that. In general, for a system with ions, both cations and anions should be included, and their charges should sum to zero.
You appear to be trying to read a restart file via IUNRE 49, but do not specify RESTART in the DYNA command.
Why read two coord sets, instead of just reading water_ions_min.crd in the first place?
Also, for a beginner, the hydronium ion seems like a poor choice; its lifetime is so short that it is not part of the standard parameter distribution.
Your HMCM command restrains the ions to a corner of the box, and not the center (0,0,0), which seems exceedingly odd.
Since you do not appear to be using the DOMDEC feature, you could use MMFP GEO PLANAR restraints to keep the ions instead the box, but the use of restraints for production MD is highly discouraged, as they perturb the system and must be accounted for during later analysis.
Overall, this is a problematic script with a problematic choice of a simulation system.
Thank you very much for your answers, Lennart and Rick! (and sorry for my late reply, I got sick)
I see I should have included more information, or removed more lines from the script. My apologies. I aimed to be brief.
I "inherited" the scrip from someone working on this before me. I changed/removed the parts I believed to already understand somewhat, one at a time, and did not dare to touch the rest yet. Hence there are fossils, like IUNRE 49 it seems, and I am thankful for any clues as to their use. I do read the output log files and search them for keywords, and I do read the docs. Nonetheless at times I can't figure out the meaning of some messages, or what to do about them. I came here because I ran out of ideas. Thank you for your valuable time!
I did read that there is no h-bond energy term in CHARMM force fields anymore, since this is taken care of by the other ff terms. The reason I use the HBOND facility is that I need to keep track of hydronium’s h-bonds - overall, I am working to investigate proton hopping. The current problem is that the HBOND facility cannot detect h-bonds across periodic boundaries. Hence the idea to keep the hydronium roughly in the box center, so its h-bonds do not escape me.
@rmv Thank you for pointing out the various details! - I will read up on the TCONS Berendsen thermostat and change it. - I am not sure I can add counter ions and change the system net charge to zero, because the ions would interact and distort h-bond patterns. - I have absolutely no idea why the coordinates are read in twice. I thought it might have something to do with CRYSTAL DEFINE…? - I simply embarrassingly misunderstood the box center position. I should have seen that, sorry. - Thanks for the suggestion of MMFP GEO PLANAR! I’ll try it out.
