Thank you Justin and Peter for your responses. I tried extending the time on the npt equilibration. It helped but not much. My final pressure after the MD run was about 1.15 bar compared to ref_p which was set to 1 bar. Peter, I will try to analyze the potential error using RMSD and Drift.
Thanks. On Mon, Apr 11, 2011 at 4:55 PM, Fabian Casteblanco < [email protected]> wrote: > Hi, > > I'm still in my first few months of using Gromacs. I started by creating > an *.itp and *.top file for *Ethanol* using CHARMM force field > parameters. I made the molecule and it looked fine, put 1000 molecules in a > box, energy minimized it to a negative potential energy, viewed it on VMD, > again looks fine. When I started running the NVT script, I set it equal to > a ref_T of 298 K. It equilibrated at the temperature. Then I tried using > an NPT script to equilibrate it to a ref_p of 1 bar. This is where I get > the problem. The output shows the density is close to the actual > experimental value of 0.789 g/cm^3. But for some reason, my pressure never > gets an average of 1 bar. It keeps oscillating, which I understand is > normal, but the average is always 1.3 or 1.4 bar (it seems the longer I let > it run, the larger the average pressure; 1.38 for 50,000 steps,dt=0.002 and > 1.45 for 75,000 steps,dt=0.002). I don't understand why the ref_p of 1 bar > is not working when I run this NPT.mdp script file. My simple goal is to > have 1000 molecules of ethanol using CHARMM ff parameters at 25degC and 1 > bar and somewhere near the experimental density. > > I would really appreciate anybody's help! I'm new to this but I'm eager to > keep getting better. > > Thanks. > > *NVT SCRIPT (this works fine and takes me to 298 K)* > File Edit Options Buffers Tools Help > title =CHARMM ETHANOL NVT equilibration > ;define =-DPOSRES ;position restrain the protein > ;Run parameters > integrator =md ;leap-frog algorithm > nsteps =50000 ;2 * 50000 = 100 ps > dt =0.002 ;2fs > ;Output control > nstxout =100 ;save coordinates every 0.2 ps > nstvout =100 ;save velocities every 0.2 ps > nstenergy =100 ;save energies every 0.2 ps > nstlog =100 ;update log file every 0.2 ps > ;Bond parameters > continuation =no ;first dynamics run > constraint_algorithm=lincs ;holonomic constraints > constraints =all-bonds ;all bonds (even heavy atom-H > bonds)constraind > lincs_iter =1 ;accuracy of LINCS > lincs_order =4 ;also related to accuracy > ;Neighborhood searching > ns_type =grid ;search neighboring grid cells > nstlist =5 ;10 fs > rlist =1.0 ;short-range neighborlist cutoff (in nm) > rcoulomb =1.0 ;short-range electrostatic cutoff (in nm) > rvdw =1.0 ;short-range van der Waals cutoff (in nm) > ;Electrostatics > coulombtype =PME ;Particle Mesh Ewald for long-range > electrostat\ > ;ics > pme_order =4 ;cubic interpolation > fourierspacing =0.16 ;grid spacing for FFT > ;Temperature coupling is on > tcoupl =V-rescale ;modified Berendsen thermostat > tc_grps =SYSTEM ;two coupling groups - more accurate > tau_t =0.1 ;0.1 ;time constant, in ps > ref_t =298 ;25 ;reference temperature, one for > each \ > ;group, in K > ;Pressure coupling is off > pcoupl =no ;no pressure coupling in NVT > ;Periodic boundary conditions > pbc =xyz ; 3-D PBC > ;Dispersion correction > DispCorr =EnerPres ;account for cut-off vdW scheme > ;Velocity generation > gen_vel =yes ;assign velocities from Maxwell > distribution > gen_temp =25 ;temperature for Maxwell distribution > gen_seed =-1 ;generate a random seed > ;END > > *NPT SCRIPT* > File Edit Options Buffers Tools Help > title =Ethanol npt equilibration > ;define =-DPOSRES ;position restrain the protein > ;Run parameters > integrator =md ;leap-frog algorithm > nsteps =50000 ;2 * 50000 = 100 ps > dt =0.002 ;2fs > ;Output control > nstxout =100 ;save coordinates every 0.2 ps > nstvout =100 ;save velocities every 0.2 ps > nstenergy =100 ;save energies every 0.2 ps > nstlog =100 ;update log file every 0.2 ps > ;Bond parameters > continuation =yes ;Restarting after NVT > constraint_algorithm=lincs ;holonomic constraints > constraints =all-bonds ;all bonds (even heavy atom-H > bonds)constraind > lincs_iter =1 ;accuracy of LINCS > lincs_order =4 ;also related to accuracy > ;Neighborhood searching > ns_type =grid ;search neighboring grid cells > nstlist =5 ;10 fs > rlist =1.0 ;short-range neighborlist cutoff (in nm) > rcoulomb =1.0 ;short-range electrostatic cutoff (in nm) > rvdw =1.0 ;short-range van der Waals cutoff (in nm) > ;Electrostatics > coulombtype =PME ;Particle Mesh Ewald for long-range > electrostat\ > ;ics > pme_order =4 ;cubic interpolation > fourierspacing =0.16 ;grid spacing for FFT > ;Temperature coupling is on > tcoupl =V-rescale ;modified Berendsen thermostat > tc-grps =SYSTEM ;two coupling groups - more accurate > tau_t =0.1; 0.1 ;time constant, in ps > ref_t =298; 300 ;reference temperature, one for > each \ > ;group, in K > ;Pressure coupling is on > pcoupl =Parrinello-Rahman ;Pressure coupling on in NPT > pcoupltype =isotropic ;uniform scaling of box vectors > tau_p =2.0 ;time constant, in ps > ref_p =1.0 ;reference pressure, in bar > compressibility =4.5e-5 ;isothermal compressibility of h2O, 1/bar > ;Periodic boundary conditions > pbc =xyz ; 3-D PBC > ;Dispersion correction > DispCorr =EnerPres ;account for cut-off vdW scheme > ;Velocity generation > gen_vel =no ;Velocity generation is off > ;gen_temp =25 ;temperature for Maxwell distribution > ;gen_seed =-1 ;generate a random seed > ;END > > > -- > *Best regards,* > ** > *Fabian F. Casteblanco* > *Rutgers University -- * > *Chemical Engineering PhD Student* > *C: +908 917 0723* > *E: **[email protected]* <[email protected]> > > -- *Best regards,* ** *Fabian F. Casteblanco* *Rutgers University -- * *Chemical Engineering PhD Student* *C: +908 917 0723* *E: **[email protected]* <[email protected]>
-- gmx-users mailing list [email protected] http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [email protected]. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
