[gmx-users] Re: energy conservation / frozen atoms
S. Alireza Bagherzadeh
s.a.bagherzadeh.h at gmail.com
Tue Jul 30 21:02:50 CEST 2013
On Tue, Jul 30, 2013 at 7:04 AM, <gmx-users-request at gromacs.org> wrote:
> Send gmx-users mailing list submissions to
> gmx-users at gromacs.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> or, via email, send a message with subject or body 'help' to
> gmx-users-request at gromacs.org
>
> You can reach the person managing the list at
> gmx-users-owner at gromacs.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of gmx-users digest..."
>
>
> Today's Topics:
>
> 1. generating user-defined topologies for surfaces (Valentina)
> 2. Re: generating user-defined topologies for surfaces
> (Justin Lemkul)
> 3. Re:
> energy conservation / frozen atoms (Justin Lemkul)
> 4. Re: generating user-defined topologies for surfaces (Valentina)
> 5. Re: generating user-defined topologies for surfaces (Valentina)
> 6. Re: Re: generating user-defined topologies for surfaces
> (Justin Lemkul)
>
>
>
>
> ------------------------------
>
> Message: 3
> Date: Tue, 30 Jul 2013 07:20:53 -0400
> From: Justin Lemkul <jalemkul at vt.edu>
> Subject: Re: [gmx-users] energy conservation / frozen atoms
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <51F7A195.80404 at vt.edu>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
>
>
> On 7/29/13 10:51 PM, S. Alireza Bagherzadeh wrote:
> > Hi All,
> >
> > I am simulating a system in which I have two solid surfaces and I keep
> them
> > frozen during simulations. I also exclude the interactions between its
> > atoms to avoid spurious contribution to the virial pressure due to large
> > forces between them as suggested in the manual.
> >
> > I run a nvt for equilibration and then I do the production run in an nve
> > ensemble; however, I am not getting good energy conservation. There is a
> > huge energy drift...
> >
> >
> > When I remove the solid surfaces, I will only have water molecules and
> > united atom methane molecules in my system. Using the same protocol I
> > obtain a very good energy conservation...
> >
>
> What happens if you unfreeze the frozen surfaces and run the same system?
>
The point is I do not want to do this as I am trying to save some
calculations.
My system is already big and including all of the bond, angles and dihedral
interactions of the solid surfaces slows down the simulation considerably.
>
> > Any insight on what might be wrong in my system would be very
> appreciated.
> >
>
> The contents of the .edr file will probably be informative, as you can
> identify
> which energy term(s) is(are) most affected. It's all probably related to
> the
> frozen surfaces themselves acting as an energy sink or something.
>
Is there anyway to fix the "energy sinking" effect of the surfaces?
I used dl_poly to simulate a similar system and there was no problem with
the energy conservation.
Many thanks,
Alireza
>
> -Justin
>
> >
> > Here is the mdp file:
> >
> > ;
> > ; File 'mdout_nve.mdp' was generated
> > ; By user: onbekend (0)
> > ; On host: onbekend
> > ; At date: Sun Jul 28 18:13:02 2013
> > ;
> >
> > ; VARIOUS PREPROCESSING OPTIONS
> > ; Preprocessor information: use cpp syntax.
> > ; e.g.: -I/home/joe/doe -I/home/mary/roe
> > include =
> > ; e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case
> sensitive)
> > define =
> >
> > ; RUN CONTROL PARAMETERS
> > integrator = md
> > ; Start time and timestep in ps
> > tinit = 0
> > dt = 0.001
> > nsteps = 1000000
> > ; For exact run continuation or redoing part of a run
> > init_step = 0
> > ; Part index is updated automatically on checkpointing (keeps files
> > separate)
> > simulation_part = 1
> > ; mode for center of mass motion removal
> > comm-mode = Linear
> > ; number of steps for center of mass motion removal
> > nstcomm = 100
> > ; group(s) for center of mass motion removal
> > comm-grps =
> >
> > ; LANGEVIN DYNAMICS OPTIONS
> > ; Friction coefficient (amu/ps) and random seed
> > bd-fric = 0
> > ld-seed = 1993
> >
> > ; ENERGY MINIMIZATION OPTIONS
> > ; Force tolerance and initial step-size
> > emtol = 10
> > emstep = 0.01
> > ; Max number of iterations in relax_shells
> > niter = 20
> > ; Step size (ps^2) for minimization of flexible constraints
> > fcstep = 0
> > ; Frequency of steepest descents steps when doing CG
> > nstcgsteep = 1000
> > nbfgscorr = 10
> >
> > ; TEST PARTICLE INSERTION OPTIONS
> > rtpi = 0.05
> >
> > ; OUTPUT CONTROL OPTIONS
> > ; Output frequency for coords (x), velocities (v) and forces (f)
> > nstxout = 0
> > nstvout = 0
> > nstfout = 0
> > ; Output frequency for energies to log file and energy file
> > nstlog = 500
> > nstcalcenergy = -1
> > nstenergy = 500
> > ; Output frequency and precision for .xtc file
> > nstxtcout = 0
> > xtc-precision = 1000
> > ; This selects the subset of atoms for the .xtc file. You can
> > ; select multiple groups. By default all atoms will be written.
> > xtc_grps = HYDW HYDG SOL GAS SiO2 SiOH
> > ; Selection of energy groups
> > energygrps = HYDW HYDG SOL GAS SiO2 SiOH
> >
> > ; NEIGHBORSEARCHING PARAMETERS
> > ; nblist update frequency
> > nstlist = 10
> > ; ns algorithm (simple or grid)
> > ns_type = grid
> > ; Periodic boundary conditions: xyz, no, xy
> > pbc = xyz
> > periodic_molecules = no
> > ; nblist cut-off
> > rlist = 1.7
> > ; long-range cut-off for switched potentials
> > rlistlong = -1
> >
> > ; OPTIONS FOR ELECTROSTATICS AND VDW
> > ; Method for doing electrostatics
> > coulombtype = PME-Switch
> > rcoulomb_switch = 1.3
> > rcoulomb = 1.5
> > ; Relative dielectric constant for the medium and the reaction field
> > epsilon_r = 1
> > epsilon_rf = 1
> > ; Method for doing Van der Waals
> > vdw-type = shift
> > ; cut-off lengths
> > rvdw-switch = 1.3
> > rvdw = 1.5
> > ; Apply long range dispersion corrections for Energy and Pressure
> > DispCorr = EnerPres
> > ; Extension of the potential lookup tables beyond the cut-off
> > table-extension = 1
> > ; Seperate tables between energy group pairs
> > energygrp_table =
> > ; Spacing for the PME/PPPM FFT grid
> > fourierspacing = 0.12
> > ; FFT grid size, when a value is 0 fourierspacing will be used
> > fourier_nx = 0
> > fourier_ny = 0
> > fourier_nz = 0
> > ; EWALD/PME/PPPM parameters
> > pme_order = 6
> > ewald_rtol = 1e-06
> > ewald_geometry = 3d
> > epsilon_surface = 0
> > optimize_fft = yes
> >
> > ; IMPLICIT SOLVENT ALGORITHM
> > implicit_solvent = No
> >
> > ; GENERALIZED BORN ELECTROSTATICS
> > ; Algorithm for calculating Born radii
> > gb_algorithm = Still
> > ; Frequency of calculating the Born radii inside rlist
> > nstgbradii = 1
> > ; Cutoff for Born radii calculation; the contribution from atoms
> > ; between rlist and rgbradii is updated every nstlist steps
> > rgbradii = 1
> > ; Dielectric coefficient of the implicit solvent
> > gb_epsilon_solvent = 80
> > ; Salt concentration in M for Generalized Born models
> > gb_saltconc = 0
> > ; Scaling factors used in the OBC GB model. Default values are OBC(II)
> > gb_obc_alpha = 1
> > gb_obc_beta = 0.8
> > gb_obc_gamma = 4.85
> > gb_dielectric_offset = 0.009
> > sa_algorithm = Ace-approximation
> > ; Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of
> GBSA
> > ; The value -1 will set default value for Still/HCT/OBC GB-models.
> > sa_surface_tension = -1
> >
> > ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> > ; Temperature coupling
> > tcoupl = no
> > nsttcouple = -1
> > nh-chain-length = 1
> > ; Groups to couple separately
> > tc_grps = system
> > ; Time constant (ps) and reference temperature (K)
> > tau_t = 0.2
> > ref_t = 370
> > ; Pressure coupling
> > Pcoupl = no
> > Pcoupltype = Isotropic
> > nstpcouple = -1
> > ; Time constant (ps), compressibility (1/bar) and reference P (bar)
> > tau_p = 0.5
> > compressibility = 4.5e-05
> > ref_p = 40.0
> > ; Scaling of reference coordinates, No, All or COM
> > refcoord_scaling = No
> > ; Random seed for Andersen thermostat
> > andersen_seed = 815131
> >
> > ; OPTIONS FOR QMMM calculations
> > QMMM = no
> > ; Groups treated Quantum Mechanically
> > QMMM-grps =
> > ; QM method
> > QMmethod =
> > ; QMMM scheme
> > QMMMscheme = normal
> > ; QM basisset
> > QMbasis =
> > ; QM charge
> > QMcharge =
> > ; QM multiplicity
> > QMmult =
> > ; Surface Hopping
> > SH =
> > ; CAS space options
> > CASorbitals =
> > CASelectrons =
> > SAon =
> > SAoff =
> > SAsteps =
> > ; Scale factor for MM charges
> > MMChargeScaleFactor = 1
> > ; Optimization of QM subsystem
> > bOPT =
> > bTS =
> >
> > ; SIMULATED ANNEALING
> > ; Type of annealing for each temperature group (no/single/periodic)
> > annealing =
> > ; Number of time points to use for specifying annealing in each group
> > annealing_npoints =
> > ; List of times at the annealing points for each group
> > annealing_time =
> > ; Temp. at each annealing point, for each group.
> > annealing_temp =
> >
> > ; GENERATE VELOCITIES FOR STARTUP RUN
> > gen_vel = no
> > gen_temp = 370
> > gen_seed = -1
> >
> > ; OPTIONS FOR BONDS
> > constraints = none
> > ; Type of constraint algorithm
> > constraint-algorithm = shake
> > ; Do not constrain the start configuration
> > continuation = no
> > ; Use successive overrelaxation to reduce the number of shake iterations
> > Shake-SOR = no
> > ; Relative tolerance of shake
> > shake-tol = 1e-10
> > ; Highest order in the expansion of the constraint coupling matrix
> > lincs-order = 4
> > ; Number of iterations in the final step of LINCS. 1 is fine for
> > ; normal simulations, but use 2 to conserve energy in NVE runs.
> > ; For energy minimization with constraints it should be 4 to 8.
> > lincs-iter = 1
> > ; Lincs will write a warning to the stderr if in one step a bond
> > ; rotates over more degrees than
> > lincs-warnangle = 30
> > ; Convert harmonic bonds to morse potentials
> > morse = no
> >
> > ; ENERGY GROUP EXCLUSIONS
> > ; Pairs of energy groups for which all non-bonded interactions are
> excluded
> > energygrp_excl = SiOH SiOH SiO2 SiO2 SiOH SiO2
> >
> > ; WALLS
> > ; Number of walls, type, atom types, densities and box-z scale factor for
> > Ewald
> > nwall = 0
> > wall_type = 9-3
> > wall_r_linpot = -1
> > wall_atomtype =
> > wall_density =
> > wall_ewald_zfac = 3
> >
> > ; COM PULLING
> > ; Pull type: no, umbrella, constraint or constant_force
> > pull = no
> >
> > ; NMR refinement stuff
> > ; Distance restraints type: No, Simple or Ensemble
> > disre = No
> > ; Force weighting of pairs in one distance restraint: Conservative or
> Equal
> > disre-weighting = Conservative
> > ; Use sqrt of the time averaged times the instantaneous violation
> > disre-mixed = no
> > disre-fc = 1000
> > disre-tau = 0
> > ; Output frequency for pair distances to energy file
> > nstdisreout = 100
> > ; Orientation restraints: No or Yes
> > orire = no
> > ; Orientation restraints force constant and tau for time averaging
> > orire-fc = 0
> > orire-tau = 0
> > orire-fitgrp =
> > ; Output frequency for trace(SD) and S to energy file
> > nstorireout = 100
> > ; Dihedral angle restraints: No or Yes
> > dihre = no
> > dihre-fc = 1000
> >
> > ; Free energy control stuff
> > free-energy = no
> > init-lambda = 0
> > delta-lambda = 0
> > foreign_lambda =
> > sc-alpha = 0
> > sc-power = 0
> > sc-sigma = 0.3
> > nstdhdl = 10
> > separate-dhdl-file = yes
> > dhdl-derivatives = yes
> > dh_hist_size = 0
> > dh_hist_spacing = 0.1
> > couple-moltype =
> > couple-lambda0 = vdw-q
> > couple-lambda1 = vdw-q
> > couple-intramol = no
> >
> > ; Non-equilibrium MD stuff
> > acc-grps =
> > accelerate =
> > freezegrps = SiO2 SiOH
> > freezedim = Y Y Y Y Y Y
> > cos-acceleration = 0
> > deform =
> >
> > ; Electric fields
> > ; Format is number of terms (int) and for all terms an amplitude (real)
> > ; and a phase angle (real)
> > E-x =
> > E-xt =
> > E-y =
> > E-yt =
> > E-z =
> > E-zt =
> >
> > ; User defined thingies
> > user1-grps =
> > user2-grps =
> > userint1 = 0
> > userint2 = 0
> > userint3 = 0
> > userint4 = 0
> > userreal1 = 0
> > userreal2 = 0
> > userreal3 = 0
> > userreal4 = 0
> >
> >
> >
> >
>
>
>
>
More information about the gromacs.org_gmx-users
mailing list