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: fabian.castebla...@gmail.com -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com
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
ethanol.itp
Description: Binary data
ethanol.top
Description: Binary data
em.gro
Description: Binary data
-- gmx-users mailing list gmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists