[gmx-users] Re: single point calculation with gromacs

[Justin Lemkul]

On 10/31/13 4:57 AM, fantasticqhl wrote:
> Dear Justin,
>
> *Thanks very much for your reply! Here is my minim.mdp I used:*
>
>
> /; minim.mdp - used as input into grompp to generate em.tpr
> ; Parameters describing what to do, when to stop and what to save
> integrator      = steep         ; Algorithm (steep = steepest descent
> minimization)
> emtol           = 1000.0        ; Stop minimization when the maximum force <
> 1000.0 kJ/mol/nm
> emstep          = 0.01          ; Energy step size
> nsteps          = 0             ; Maximum number of (minimization) steps to
> perform
> ; Parameters describing how to find the neighbors of each atom and how to
> calculate the interactions
> nstlist         = 1             ; Frequency to update the neighbor list and
> long range forces
> ns_type         = simple                ; Method to determine neighbor list
> (simple, grid)
> rlist           = 9.0           ; Cut-off for making neighbor list (short
> range forces)
> coulombtype     = Cut-off               ; Treatment of long range
> electrostatic interactions
> rcoulomb        = 9.0           ; Short-range electrostatic cut-off
> rvdw            = 9.0           ; Short-range Van der Waals cut-off
> pbc             = no            ; Periodic Boundary Conditions (yes/no)
> /
>
>
> *If I used this mdp file for 0-step minimization, I could get potential
> energy around 700 kJ/mol for my system. But I would get potential energy
> around 2.6e+05 kJ/mol if I used the mdv.mdp (md simulation in vacuum) for
> calculations with rerun option. And the bond potential could also reach as
> high as 1.1e+05 kJ/mol. Actually, my system is very small, it contains only
> 37 atoms. I believe that the energy reported by 0-step minimization were
> more reasonable. So I guess that there might be some problem for the mdp
> file usd with rerun, and here is the mdv.mdp:*
>
>
>
> /define              = ;-DPOSRES
> integrator          =  md       ; molecular dynamics algorithm
> tinit               =  0.0      ; start time and timestep in ps
> dt                  =  0.002    ; time step in ps
> nsteps              =  2        ; number of steps for 1000ns run
> emtol               =  100    ; convergence criterion
> emstep              =  0.05      ; intial step size
> nstlist             =  0       ; step frequency for updating neighbour list
> ns_type             =  simple     ; method for neighbour searching (?)
> nstxout             =  1    ; frequency for writing coords to output .trr
> file
> nstvout             =  1     ; frequency for writing velocities to
> output...should be same as nstxout
> nstfout             =  1        ; frequency for writing forces to output
> nstlog              =  1      ; frequency for writing energies to log file
> nstenergy           =  1      ; frequency for writing energies to energy
> file
> nstxtcout           =  1     ; frequency for writing coords to xtc traj
> xtc_grps            =  system   ; group(s) whose coords are to be written in
> xtc traj
> energygrps          =  system   ; group(s) whose energy is to be written in
> energy file
> pbc                 =  no      ; use pbc
> rlist               =  0      ; cutoff lengths (nm)
> epsilon_r           =  1.0      ; Dielectric constant (DC) for twin-range or
> DC of reaction field
> niter               =  100      ; Some thingies for future use
> fourierspacing      =  0.16
> fourier_nx          =  30
> fourier_ny          =  30
> fourier_nz          =  30
> coulombtype         =  Cut-off      ; truncation for minimisation, with
> large cutoff
> rcoulomb            =  0
> rcoulomb-switch     =  0
> vdw-type                 = Cut-off  ; truncation for minimisation, with
> large cutoff
> rvdw-switch              = 0
> rvdw                     = 0   ; cut-off lengths
> epsilon_surface          = 0
> optimize_fft             = yes
> Tcoupl              =  V-rescale
> tc_grps             = system
> tau_t               = 0.01
> ref_t               = 300
> Pcoupl              = no ; Parrinello-Rahman ; Pressure coupling
> gen_vel             =  yes
> gen_temp            =  300
> gen_seed            =  -1
> constraints         = none  ; OPTIONS FOR BOND CONSTRAINTS
> constraint-algorithm  = Lincs   ; Type of constraint algorithm
> lincs_order         =  4        ;4        ; Highest order in the expansion
> of the constraint coupling matrix
> lincs_iter          =  1
> lincs_warnangle     =  30       ; Lincs will write a warning to the stderr
> if in one step a bond rotates
>                                  ; over more degrees than
> unconstrained-start      = no   ; Do not constrain the start configuration
> ;Shake-SOR                = no   ; Use successive overrelaxation to reduce
> the number of shake iterations
> ;shake-tol                = 1e-04 ; Relative tolerance of shake
> morse                    = no   ; Convert harmonic bonds to morse potentials
> /
>
>
> *Could you please have a check for me again? Thanks in advance!
>

Was anything in the bug report I posted before relevant?  If you're using any 
version in the 4.6.x series, you may be observing some of the same behavior I did.

There are numerous differences between these .mdp files.  You should stay 
consistent in terms of cutoff length, application of constraints, etc.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441

==================================================