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

Francesca Vitalini francesca.vitalini11 at gmail.com
Thu Jan 23 10:14:37 CET 2014


Hi,

thank you for the response,


Il giorno 23/gen/2014, alle ore 06:35, bipin singh <bipinelmat at gmail.com> ha scritto:

>> 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.
> 

If I transform a HG into a HB then won't I have problems with the bonds definitions? How should I change my topology in order to define that in the end state the mutated HG has to be bound to CB instead?

Thank you very much.

Francesca

> 
> 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*
> -- 
> 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.



More information about the gromacs.org_gmx-users mailing list