[gmx-users] NPT for compressing the system

[Mark Abraham]

On 14/07/2010 1:09 PM, Moeed wrote:
> Hello,
>
> Thank you for suggestions. Yes simulation crashes. I am not able see 
> the whole run (1000ps). For all pressures I tired the final simulation 
> box is 2.2*2.2*2.2 nm but simulation crashes before 1000ps. The 
> longest simulation is for 120 bar (727ps). and I think that is why I 
> am not getting the _after_md.gro file after mdrun in none of the 
> simulations. I wanted to take that convoluted structure and replicate 
> the box.. 

I don't recall why you have a large box that you have to compress, but 
the way to do that is to be gentler with it. What you're doing is 
integrating equations of motion with forces that are much larger than 
those normally experienced, and that is not a numerically stable 
situation if you keep the timestep the same size, since the atoms will 
travel so far that the approximations in the integration scheme are no 
longer reasonable. Think about how (badly) a similar model of a simple 
pendulum would work if the timestep was comparable with the period of 
the motion, and someone was bashing the pendulum each time step...

A few bars (maybe 10?) of overpressure will compress the box, and do so 
in such a way that the system does not get too far from equilbrium. That 
maximizes your chance of reaching your desired box size with a system 
that is not greatly unhappy. Even if such an unhappy system was not 
exploding, you'll still have to equilibrate such a system for longer 
than one that had a more gentle compression regime. Write the 
coordinates fairly often, and then even if you overshoot the 
compression, you can get a structure of about the right size to start 
again to equilibrate at about the right size with 1 bar P-coupling.

There's no value in doing the compression fast if fast doesn't work!

Also, generating velocities and then immediately trying to compress the 
box before anything has equilibrated is asking for trouble. If your 
system density is very far from equilibrium, you might need to use a 
small timestep to get anything to work

Mark

> mdrun -s PE60-P120_md.tpr -o PE60-P120_md.tpr -c PE60-P120_after_md -v>&  output.mdrun_md
> As you suggested I do the runs in stepwise manner. I need still a 
> little smaller box size than 2*2*2 but I noticed sth during these 
> simulations. When I view the trajectory in ngmx on the menu bar I can 
> select type of box. When I select triclinic or octahedron the final 
> confromation of the chain is more convoluted and fits completely in 
> the box while in case of rectangular box except for pressure of 120 
> bar for all other pressures I tried most parts of the chain remain 
> outside of the small box. molecules is more extended...I am wondeirng 
> if positions of atoms are calculated and stored how changing the box 
> type in trajectory viewer can affect the structure?!
> ***********************************************************************************************mdp 
> file I used:
> title               =  PE-Hexane
> ;define           =  -DPOSRES    ; tells gromacs to perform position 
> restrained dynamics/include posre.itp into topology used for position 
> restraint
> pbc              =  xyz        ; use priodic BCs in all directions
>
> ;        Run control
> integrator          =  md        ; type of dynamics algorithm. Here md 
> uses a leap-frog algorithm for integrating Newtons's eq of motion
> dt                  =  0.002        ; in ps !
> nsteps              =  500000        ; length of simulation= nsteps*dt
> nstcomm             =  1        ; frequency for center of mass motion 
> removal
>
> ;        Output control
> nstenergy           =  100        ; frequency to write energies to 
> energy file. i.e., energies and other statistical data are stored 
> every 10 steps
> nstxout             =  100        ; frequency to write 
> coordinates/velocity/force to output trajectory file. how often 
> snapshots are collected= nstxout*dt
> nstvout             =  100
> nstfout             =  0
> nstlog              =  100        ; frequency to write energies to log 
> file
> nstxtcout          =  10        ; frequency to write coordinates to 
> xtc trajectory
>
> ;        Neighbor searching
> nstlist             =  10        ; frequency to update neighbor list. 
> Neighborlist will be updated at least every 10 steps. Manual p80
> ns_type             =  grid        ; make a grid in the box and only 
> check atoms in neighboring grid cells when constructing a new neighbor 
> list every nstlist steps
>
> ;        Electrostatics/VdW
> coulombtype         =  Shift        ; tells gromacs how to model 
> electrostatics. Coulomb/LJ potential is decreased over the whole range 
> and forces decay smoothly to zero between
> vdw-type            =  Shift      ; rcoulomb-switch/rvw-switch & 
> rcoulomb/rvdw
> rcoulomb-switch     =  0        ; where to start switching the Coulomb 
> potential
> rvdw-switch         =  0        ; where to start switching the LJ 
> potential
>
> ;        Cut-offs
> rlist               =  1.1        ; in nm. Cut-off distance for 
> short-range neighbor list
> rcoulomb            =  1.0        ; distance for coulomb cut-off
> rvdw                =  1.0        ; distance for coulomb cut-off
>
> ;        Temperature coupling
> Tcoupl              =  berendsen
> tc-grps             =  System  ;HEX            ; groups to couple to 
> thermostat; Berendsen temperature coupling is on in these groups
> tau_t               =  0.1     ;0.1            ; time constant for T 
> coupling
> ref_t               =  300     ;300            ; reference T for 
> coupling. When you alter the T, don't forget to change the gen_temp 
> for velocity generation
>
> ;        Pressure coupling
> Pcoupl              =  berendsen    ; Pressure coupling is not on
> Pcoupltype          =  isotropic     ; means the box expands and 
> contracts in all directions (x,y,z) in order to maintain the proper 
> pressure
> tau_p               =  0.5        ; time constant for coupling in ps
> compressibility     =  4.5e-5        ; compressibility of solvent used 
> in simulation in 1/bar
> ref_p               =  160.0        ; reference P for coupling in bar
>
> ;        Velocity generation      Generate velocites is on at 300 K. 
> Manual p155
> gen_vel             =  yes        ; generate velocites according to 
> Maxwell distribution at T: gen_temp with random gen seed gen_seed
> gen_temp            =  300.0        ; T for Maxwell distribution
> gen_seed            =  173529        ; used to initialize random 
> generator for random velocities
>
> ;        Bonds
> constraints         =  all-bonds    ; sets the LINCS constraint for 
> all bonds
> constraint-algorithm = lincs
> *******************************************************************************************
>
>
> Thanks
>
>
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100714/c013cba0/attachment.html>