[gmx-users] FEP

Fabian Casteblanco fabian.casteblanco at gmail.com
Thu Apr 26 17:43:23 CEST 2012


Hello all,

This is in reply to Michael shirts a while ago on a FEP of a R-CH3 to
an R-H group.  Below is the orignal email.  I recently tested out
mutating a CH3-CH3-(3dummy atoms) molecule on both sides in order to
test out that a peturbation would give you a total of 0.  forcefield
used was CGenFF

So for procedure was that I turned the CH3 group on the left to an H
with 3 dummy atoms and I turned an H and 3 dummy atoms on the right to
a CH3 group, essentially mutating the molecule to another version of
itself (should give you a total dG of 0 if it works).

STEP 1. (uncharging everything except the -CH2 group inbetween both
sides  ....   CH3-CH2-HD3)  *D are dummy atoms
 [ atoms ]
;  nr  type  resnr  resid    atom  cgnr  charge    mass    total_charge
    1    CG331     1   SIM       C1     1   	 -0.27
	12.0110		CG331	0.00	 12.0110
    2     HGA3     1   SIM      H11     2    	  0.09
	1.00800		HGA3	0.00	 1.00800
    3     HGA3     1   SIM      H12     3    	  0.09
	1.00800		HGA3	0.00	 1.00800
    4     HGA3     1   SIM      H13     4    	  0.09
	1.00800		HGA3	0.00	 1.00800
    5    CG331     1   SIM       C2     5	 	 -0.27 	12.0110		CG331
-0.18	 12.0110
    6	  HGA3	 1   SIM      H21     6	  	  0.09   			
    7	  HGA3	 1   SIM      H22     7    	  0.09 				
    8	  HGA3	 1   SIM      H23     8    	  0.09   	1.00800		HGA3	0.00	 1.00800		
    9	   DUM	 1   SIM      H31     9	  	  0.00   	1.00800						
   10	   DUM	 1   SIM      H32    10     	  0.00 	1.00800						
   11	   DUM	 1   SIM      H33    11    	  0.00   	1.00800		

STEP 2.  (mutation of groups.   	CH3-CH2-HD3 becomes D3H-CH2-CH3)
[ atoms ]
;  nr  type  resnr  resid    atom  cgnr  charge    mass    total_charge
    1    CG331      1   SIM       C1      1   	  0.00
	12.0110		HGA3	0.00	   1.00800
    2     HGA3      1   SIM      H11      2    	  0.00
	1.00800		DUM	        0.00	   1.00800
    3     HGA3      1   SIM      H12      3    	  0.00
	1.00800		DUM	        0.00	   1.00800
    4     HGA3      1   SIM      H13      4    	  0.00  	1.00800		DUM	
       0.00	   1.00800
    5    CG331      1   SIM       C2      5	 	 -0.18 	  	
    6	  HGA3	  1   SIM      H21     6	  	  0.09   	
    7	  HGA3	  1   SIM      H22     7    	  0.09 		
    8	  HGA3	  1   SIM      H23     8    	  0.00
	1.00800		CG331	0.00	   12.0110
    9	   DUM	  1   SIM      H31     9	  	  0.00   	1.00800		HGA3	0.00	
  1.00800
   10	   DUM	  1   SIM      H32    10     	  0.00 	1.00800		HGA3	0.00	
  1.00800
   11	   DUM	  1   SIM      H33    11    	  0.00
	1.00800		HGA3	0.00	   1.00800

STEP3.  (opposite of step1)
[ atoms ]
;  nr  type  resnr  resid    atom  cgnr  charge    mass    total_charge
    1     HGA3   1   SIM         C1     1   	  0.00 		12.0110		HGA3
  0.09	  1.00800
    2      DUM   1   SIM        H11     2    	  0.00   	
1.00800		DUM	       0.00	  1.00800
    3      DUM   1   SIM        H12     3    	  0.00   	
1.00800		DUM	       0.00	  1.00800
    4      DUM   1   SIM        H13     4    	  0.00  	
1.00800		DUM	       0.00	  1.00800
    5    CG331   1   SIM         C2     5	 	 -0.18 	  	12.0110		CG331
  -0.27	  12.0110
    6	  HGA3	 1   SIM      H21     6	  	  0.09   	
    7	  HGA3	 1   SIM      H22     7    	  0.09 		
    8	 CG331	 1   SIM      H23     8    	  0.00   	
1.00800		CG331     -0.27	  12.0110
    9	  HGA3	 1   SIM      H31     9	  	  0.00   	
1.00800		HGA3	0.09	  1.00800
   10	  HGA3	 1   SIM      H32    10     	  0.00 		
1.00800		HGA3	0.09	  1.00800
   11	  HGA3	 1   SIM      H33    11    	  0.00   	
1.00800		HGA3	0.09	  1.00800

So in the end, the procedure worked using the g_bar method.  I took 20
simulations for steps 1 and 3 and 50 simulations for step 2.  Here are
the results:

Step1: -37.62 +/- 0
Step2:  -0.12 +/- 0.24
Step3:  37.68 +/- 0.13
TOTAL ~ -0.06 kJ/mol  which is expected

My question is that I did something similar for a larger molecule (65
atoms) and someone said before that step 1 and 3 should be very small.
 I'm getting values in the near -20 kJ/mol for step 1 and 3.  I feel
it should be straightforward since I'm simply uncharging a CH3 group
and I don't get why I'm getting such a large number.  If it is suppose
to be very small, is it safe to assume that step2 should be a majority
of the final value.

If anyone has any experience in this area all help is greatly
appreciated.  Thank you very much.

-Fabian




> So BAR is only
> referring to the mathematical code used to calculate the overall free
> energy for the FEP, correct?

Yes.  The information required is the same, with the exception that
exponential averaging requires energy differences from only one state,
whereas BAR requires energy differences from both directions.  Since
you are likely running simulations at a range of states anyway,
there's no extra cost.

> My question is, for a mutation of a -CH3 group to a -H group, is it
> better to simply run:
> [+ from (Lambda=0 ,  R-CH3, full charges and interactions -STATE A)
> --> (Lambda=1, R-CH, full charges and interactions -STATE B)]
>
> OR
>
> [1) from (Lambda=0 ,  R-CH3, STATE A : Charges and LJ Interactions: ON)
> 2)  (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated)
> 3)  (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF)
> 4)  (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to -H)
> 5)  (Lambda=1, R-CH3, -CH3  STATE B : Charges and LJ Interactions: ON)

Currently, you'd need to run three simulations to do this (will likely
be simplified somewhat in 4.6).

Set 1: turn R-CH3 charges off in a way that preserves the total charge.
Set 2: change CH3 LJ to H LJ
Set 3: Turn R-H charges on in a way that preserves the total charge.

Note that the first and last transformations will need only one two
states -- their energies will be VERY small.

> Reason I'm asking is because when I try the first choice to do it
> STATE A to STATE B in one step, when I reach Lambda=0.85 and above on
> the NVT equilibration right after EM, I receive errors saying that
> bonds are moving way to far off their constraints which leads me to
> believe that the system is moving too far from where it was energy
> minimized.  Errors such as:
>
> Step 188, time 0.376 (ps)
> LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.000017, max 0.000636 (between atoms 9 and 68)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>      9     68   31.2    0.1111   0.1110      0.1111
>
> Step 188, time 0.376 (ps)  LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.000015, max 0.000531 (between atoms 9 and 68)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>      9     68   31.0    0.1111   0.1110      0.1111

You probably have problems with the charge + LJ terms on the hydrogen.
 Note that you could possibly use soft core for the Coloumb
interactions for hydrogen -- it's sometimes causes sampling problems,
but for H's it probably doesn't matter.  I haven't done this approach
much so, so I can't tell for sure if it will work.


--
Best regards,

Fabian F. Casteblanco
Rutgers University --
Chemical Engineering



More information about the gromacs.org_gmx-users mailing list