[gmx-users] regarding free energy perturbation

[zhangw]

Dear Gromacs users,
 I am using Gromacs 5.0.6 with CHARMM27 forcefield. I want to do FEP analysis and I did it based on this tutorials: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/01_theory.html

. 
But I just want to peturb some atoms of the protein, however, in the mdp file,I only can choose whole protein moleculetype.   couple-moltype = Protein
I don’t know how to calculate free energy by only perturbing charge of some atoms of protein. 
and I have tried to modify topology file  that I added the state B, liking the following: but it seems unable to recognize the state B in topology file.

    25          C      3    ALA      C     25       0.51     12.011  C  0.00  12.011 ; qtot 0.465
    26          O      3    ALA      O     26      -0.51     15.999  O  0.00  15.999 ; qtot -0.045
; residue   4 ALA rtp ALA  q  0.0
    27        NH1      4    ALA      N     27      -0.47     14.007  NH1 0.00  14.007 ; qtot -0.515
    28          H      4    ALA     HN     28       0.31      1.008  H  0.00  1.008 ; qtot -0.205
    29        CT1      4    ALA     CA     29       0.07     12.011  CT1 0.00  12.011 ; qtot -0.135
    30         HB      4    ALA     HA     30       0.09      1.008  HB  0.00  1.008 ; qtot -0.045
    31        CT3      4    ALA     CB     31      -0.27     12.011   ; qtot -0.315
    32         HA      4    ALA    HB1     32       0.09      1.008   ; qtot -0.225
    33         HA      4    ALA    HB2     33       0.09      1.008   ; qtot -0.135
    34         HA      4    ALA    HB3     34       0.09      1.008   ; qtot -0.045
    35          C      4    ALA      C     35       0.51     12.011   ; qtot 0.465
    36          O      4    ALA      O     36      -0.51     15.999   ; qtot -0.045
; residue   5 ALA rtp ALA  q  0.0
    37        NH1      5    ALA      N     37      -0.47     14.007   ; qtot -0.515
    38          H      5    ALA     HN     38       0.31      1.008   ; qtot -0.205


but it seems that it can't recognize state B in topology file. And it seems  to calculate the free energy just based on the mdp file, couple-moltype = Protein
Mdp file is the following:

; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 1.0
compressibility          = 4.5e-05
ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda_state        = 0
delta_lambda             = 0
calc_lambda_neighbors    = 1        ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state        0    1    2    3    4    5    6    7    8    9    10   11
vdw_lambdas              = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
coul_lambdas             = 0.0 0.01 0.02 0.04 0.06 0.08 0.1 0.15 0.25 0.5 0.75 1.0
; We are not transforming any bonded or restrained interactions
bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha                 = 0.5
sc-coul                  = no       ; linear interpolation of Coulomb (none in this case)
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-moltype           = Protein  ; name of moleculetype to decouple
couple-lambda0           = vdw-q     ; only van der Waals interactions
couple-lambda1           = none     ; turn off everything, in this case only vdW
couple-intramol          = no
nstdhdl                  = 10
; Do not generate velocities
gen_vel                  = no
; options for bonds
constraints              = h-bonds  ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm     = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation             = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 12


I’m wandering how can I solve the problem ? It would be a great help if anyone could help me out here.
Many thanks!

Zhang Xiao