[gmx-users] Free energy perturbation parameters sc-alpha sc-power sc-sigma

bipin singh bipinelmat at gmail.com
Thu Jan 23 06:35:42 CET 2014


>I am trying to perform a series of simulations where I gradually transform
>a valine capped amino acid (Ac-VAL-NHMe) into an alanine capped amino acid
>(Ac-ALA-NHMe), ie I have to turn on the HG and transform CG into HB. I have
>defined 6 lambda values at which I would like to equilibrate my system
>(lambda=0,0.2,0.40.6,0.8,1)

Would it not be better if you turn off the vdw and coul. terms of  CG (1
atom) and HG (two atoms) and transform one of the HG to HB while doing the
transformation from valine to alanine, rather than mutating the CG to HB.


On Wed, Jan 22, 2014 at 9:40 PM, Francesca Vitalini <
francesca.vitalini11 at gmail.com> wrote:

> Dear Gromacs users,
>
> I am trying to perform a series of simulations where I gradually transform
> a valine capped amino acid (Ac-VAL-NHMe) into an alanine capped amino acid
> (Ac-ALA-NHMe), ie I have to turn on the HG and transform CG into HB. I have
> defined 6 lambda values at which I would like to equilibrate my system
> (lambda=0,0.2,0.40.6,0.8,1)
>
>  I am using GROMACS.4.5.5. and OPLS-AA force field.
>
> I have been following the tutorial by Justin Lemkul
>
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html
>
> I have problems with the nvt equilibration with lambda=0.8 and lambda=1 as
> my systems systematically crashes. I believe that the reason is in the
> tuning of the parameters sc-alpha sc-power sc-sigma in the mdp file.
> I understand from the manual that sc_sigma should not be the default option
> as I have 6 hydrogens disappearing, and consequently the sc-alpha sc-power
> parameters will be effected as the interactions of the hydrogens in the
> soft core state will get stronger.
> However I have no experience in how to select these parameters and would
> like to ask, if anyone has more experience, which could be plausible
> values.
>
> Thank you very much,
>
> Francesca
>
>
> For completeness I attach here a copy of my topology and mdp files.
>
>
> TOPOLOGY
>
>
> ; Include forcefield parameters
> #include "oplsaa.ff/forcefield.itp"
>
> [ moleculetype ]
> ; Name            nrexcl
> Protein_chain_A     3
>
> [ atoms ]
> ;   nr       type  resnr residue  atom   cgnr     charge       mass
> typeB    chargeB      massB
> ; residue   1 ACE rtp ACE  q  0.0
>      1   opls_135      1    ACE    CH3      1      -0.18     12.011
> opls_135     -0.18     12.011
>      2   opls_140      1    ACE   HH31      1       0.06      1.008
> opls_140      0.06      1.008
>      3   opls_140      1    ACE   HH32      1       0.06      1.008
> opls_140      0.06      1.008
>      4   opls_140      1    ACE   HH33      1       0.06      1.008
> opls_140      0.06      1.008
>      5   opls_235      1    ACE      C      2        0.5     12.011
> opls_235       0.5     12.011
>      6   opls_236      1    ACE      O      2       -0.5    15.9994
> opls_236      -0.5    15.9994
> ; residue   2 VAL rtp VAL  q  0.0
>      7   opls_238      2    VAL      N      3       -0.5    14.0067
> opls_238      -0.5    14.0067
>      8   opls_241      2    VAL      H      3        0.3      1.008
> opls_241       0.3      1.008
>      9  opls_224B      2    VAL     CA      3       0.14     12.011
> opls_224B      0.14     12.011
>     10   opls_140      2    VAL     HA      3       0.06      1.008
> opls_140      0.06      1.008
>     11   opls_137      2    VAL     CB      4      -0.06     12.011
> opls_137     -0.18     12.011
>     12   opls_140      2    VAL     HB      4       0.06      1.008
> opls_140      0.06      1.008
>     13   opls_135      2    VAL    CG1      5      -0.18     12.011
> opls_140      0.06      1.008
>     14   opls_140      2    VAL   HG11      5       0.06      1.008
> opls_140      0.00      0.000
>     15   opls_140      2    VAL   HG12      5       0.06      1.008
> opls_140      0.00      0.000
>     16   opls_140      2    VAL   HG13      5       0.06      1.008
> opls_140      0.00      0.000
>     17   opls_135      2    VAL    CG2      6      -0.18     12.011
> opls_140      0.06      1.008
>     18   opls_140      2    VAL   HG21      6       0.06      1.008
> opls_140      0.00      0.000
>     19   opls_140      2    VAL   HG22      6       0.06      1.008
> opls_140      0.00      0.000
>     20   opls_140      2    VAL   HG23      6       0.06      1.008
> opls_140      0.00      0.000
>     21   opls_235      2    VAL      C      7        0.5     12.011
> opls_235       0.5     12.011
>     22   opls_236      2    VAL      O      7       -0.5    15.9994
> opls_236      -0.5    15.9994
> ; residue   3 NAC rtp NAC  q  0.0
>     23   opls_238      3    NAC      N      8       -0.5    14.0067
> opls_238      -0.5    14.0067
>     24   opls_241      3    NAC      H      8        0.3      1.008
> opls_241       0.3      1.008
>     25   opls_242      3    NAC    CH3      9       0.02     12.011
> opls_242      0.02     12.011
>     26   opls_140      3    NAC   HH31      9       0.06      1.008
> opls_140      0.06      1.008
>     27   opls_140      3    NAC   HH32      9       0.06      1.008
> opls_140      0.06      1.008
>     28   opls_140      3    NAC   HH33      9       0.06      1.008
> opls_140      0.06      1.008
>
> [ bonds ]
> ;  ai    aj funct            c0            c1            c2            c3
>     1     2     1
>     1     3     1
>     1     4     1
>     1     5     1
>     5     6     1
>     5     7     1
>     7     8     1
>     7     9     1
>     9    10     1
>     9    11     1
>     9    21     1
>    11    12     1
>    11    13     1
>    11    17     1
>    13    14     1
>    13    15     1
>    13    16     1
>    17    18     1
>    17    19     1
>    17    20     1
>    21    22     1
>    21    23     1
>    23    24     1
>    23    25     1
>    25    26     1
>    25    27     1
>    25    28     1
>
> [ pairs ]
> ;  ai    aj funct            c0            c1            c2            c3
>     1     8     1
>     1     9     1
>     2     6     1
>     2     7     1
>     3     6     1
>     3     7     1
>     4     6     1
>     4     7     1
>     5    10     1
>     5    11     1
>     5    21     1
>     6     8     1
>     6     9     1
>     7    12     1
>     7    13     1
>     7    17     1
>     7    22     1
>     7    23     1
>     8    10     1
>     8    11     1
>     8    21     1
>     9    14     1
>     9    15     1
>     9    16     1
>     9    18     1
>     9    19     1
>     9    20     1
>     9    24     1
>     9    25     1
>    10    12     1
>    10    13     1
>    10    17     1
>    10    22     1
>    10    23     1
>    11    22     1
>    11    23     1
>    12    14     1
>    12    15     1
>    12    16     1
>    12    18     1
>    12    19     1
>    12    20     1
>    12    21     1
>    13    18     1
>    13    19     1
>    13    20     1
>    13    21     1
>    14    17     1
>    15    17     1
>    16    17     1
>    17    21     1
>    21    26     1
>    21    27     1
>    21    28     1
>    22    24     1
>    22    25     1
>    24    26     1
>    24    27     1
>    24    28     1
>
> [ angles ]
> ;  ai    aj    ak funct            c0            c1
> c2            c3
>     2     1     3     1
>     2     1     4     1
>     2     1     5     1
>     3     1     4     1
>     3     1     5     1
>     4     1     5     1
>     1     5     6     1
>     1     5     7     1
>     6     5     7     1
>     5     7     8     1
>     5     7     9     1
>     8     7     9     1
>     7     9    10     1
>     7     9    11     1
>     7     9    21     1
>    10     9    11     1
>    10     9    21     1
>    11     9    21     1
>     9    11    12     1
>     9    11    13     1
>     9    11    17     1
>    12    11    13     1
>    12    11    17     1
>    13    11    17     1
>    11    13    14     1
>    11    13    15     1
>    11    13    16     1
>    14    13    15     1
>    14    13    16     1
>    15    13    16     1
>    11    17    18     1
>    11    17    19     1
>    11    17    20     1
>    18    17    19     1
>    18    17    20     1
>    19    17    20     1
>     9    21    22     1
>     9    21    23     1
>    22    21    23     1
>    21    23    24     1
>    21    23    25     1
>    24    23    25     1
>    23    25    26     1
>    23    25    27     1
>    23    25    28     1
>    26    25    27     1
>    26    25    28     1
>    27    25    28     1
>
> [ dihedrals ]
> ;  ai    aj    ak    al funct            c0            c1
> c2            c3            c4            c5
>     2     1     5     6     3
>     2     1     5     7     3
>     3     1     5     6     3
>     3     1     5     7     3
>     4     1     5     6     3
>     4     1     5     7     3
>     1     5     7     8     3
>     1     5     7     9     3
>     6     5     7     8     3
>     6     5     7     9     3
>     5     7     9    10     3
>     5     7     9    11     3
>     5     7     9    21     3
>     8     7     9    10     3
>     8     7     9    11     3
>     8     7     9    21     3
>     7     9    11    13     3    dih_VAL_chi1_N_C_C_C
>     7     9    11    17     3    dih_VAL_chi1_N_C_C_C
>    21     9    11    13     3    dih_VAL_chi1_C_C_C_CO
>    21     9    11    17     3    dih_VAL_chi1_C_C_C_CO
>     7     9    11    12     3
>    10     9    11    12     3
>    10     9    11    13     3
>    10     9    11    17     3
>    21     9    11    12     3
>     7     9    21    22     3
>     7     9    21    23     3
>    10     9    21    22     3
>    10     9    21    23     3
>    11     9    21    22     3
>    11     9    21    23     3
>     9    11    13    14     3
>     9    11    13    15     3
>     9    11    13    16     3
>    12    11    13    14     3
>    12    11    13    15     3
>    12    11    13    16     3
>    17    11    13    14     3
>    17    11    13    15     3
>    17    11    13    16     3
>     9    11    17    18     3
>     9    11    17    19     3
>     9    11    17    20     3
>    12    11    17    18     3
>    12    11    17    19     3
>    12    11    17    20     3
>    13    11    17    18     3
>    13    11    17    19     3
>    13    11    17    20     3
>     9    21    23    24     3
>     9    21    23    25     3
>    22    21    23    24     3
>    22    21    23    25     3
>    21    23    25    26     3
>    21    23    25    27     3
>    21    23    25    28     3
>    24    23    25    26     3
>    24    23    25    27     3
>    24    23    25    28     3
>
> [ dihedrals ]
> ;  ai    aj    ak    al funct            c0            c1
> c2            c3
>     1     7     5     6     1    improper_O_C_X_Y
>     5     9     7     8     1    improper_Z_N_X_Y
>     9    23    21    22     1    improper_O_C_X_Y
>    21    25    23    24     1    improper_Z_N_X_Y
>
> ; Include Position restraint file
> #ifdef POSRES
> #include "Ac_V_NHMe/Ac_V_NHMe.itp"
> #endif
>
> ; Include water topology
> #include "oplsaa.ff/tip3p.itp"
>
> #ifdef POSRES_WATER
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ;  i funct       fcx        fcy        fcz
>    1    1       1000       1000       1000
> #endif
>
> ; Include topology for ions
> #include "oplsaa.ff/ions.itp"
>
> [ system ]
> ; Name
> ACETYL-VALINE-METHYLAMIDE in water
>
> [ molecules ]
> ; Compound        #mols
> Protein_chain_A     1
> SOL               682
>
>
>
>
> NVT_VDW-Q.mdp
> here I intend to switch off the coulomb interactions.
>
> ; Run control
> integrator               = sd       ; Langevin dynamics
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 10000    ; 20 ps
> nstcomm                  = 100
> ; Output control
> nstxout                  = 10000
> nstvout                  = 10000
> nstfout                  = 0
> nstlog                   = 1000
> nstenergy                = 500
> nstxtcout                = 500
> xtc-precision            = 500
> xtc_grps                 = Protein
> ; Neighborsearching and short-range nonbonded interactions
> nstlist                  = 5
> ns_type                  = grid
> pbc                      = xyz
> rlist                    = 1.0
> ; Electrostatics
> coulombtype              = PME
> rcoulomb                 = 1.0
> ; van der Waals
> vdw-type                 = switch
> rvdw-switch              = 0.8
> rvdw                     = 0.9
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr                  = EnerPres
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing           = 0.12
> ; EWALD/PME/PPPM parameters
> pme_order                = 6
> ewald_rtol               = 1e-06
> epsilon_surface          = 0
> optimize_fft             = no
> ; Temperature coupling
> ; tcoupl is implicitly handled by the sd integrator
> tc_grps                  = system
> tau_t                    = 1.0
> ref_t                    = 300
> ; Pressure coupling is off for NVT
> Pcoupl                   = No
> tau_p                    = 0.5
> compressibility          = 4.5e-05
> ref_p                    = 1.0
> ; Free energy control stuff
> free_energy              = yes
> init_lambda              = 0.8
> delta_lambda             = 0
> sc-alpha                 = 0.5
> sc-power                 = 1.0
> sc-sigma                 = 0.005
> couple-moltype           = Protein_chain_A      ; name of moleculetype to
> decouple
> couple-lambda0           = vdw-q      ; only van der Waals interactions
> couple-lambda1           = vdw     ; turn off everything, in this case only
> vdW
> couple-intramol          = no
> ; Generate velocities to start
> gen_vel                  = yes
> gen_temp                 = 300
> gen_seed                 = -1
> ; options for bonds
> constraints              = h-bonds  ; we only have C-H bonds here
> ; Type of constraint algorithm
> constraint-algorithm     = lincs
> ; Do not constrain the starting configuration
> continuation             = no
> ; Highest order in the expansion of the constraint coupling matrix
> lincs-order              = 12
>
>
> NVT_VDW.mdp
> here I intend to switch off the Van der Waals interactions. The only
> difference with the previous file is in the Free energy parameters, so I
> will report only those.
>
> ; Free energy control stuff
> free_energy              = yes
> init_lambda              = 0.8
> delta_lambda             = 0
> sc-alpha                 = 0.5
> sc-power                 = 1.0
> sc-sigma                 = 0.3
> couple-moltype           = Protein_chain_A      ; name of moleculetype to
> decouple
> couple-lambda0           = vdw      ; only van der Waals interactions
> couple-lambda1           = none     ; turn off everything, in this case
> only vdW
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>



-- 



*--------------------Thanks and Regards,Bipin Singh*


More information about the gromacs.org_gmx-users mailing list