[gmx-users] Re: FEP

[Fabian Casteblanco]

Hello Michael or others,

It seems I'm still getting errors when doing a FEP on a molecule (a
-CH3 to a -H).  This below is for when I was charging things from a -H
uncharged to -H charged, although it also happens when I'm actually
converting the -CH3 to -H (at Lambdas greater than 85%).  I made sure
to keep the charges balanced at 0 while mutating and I did it at 3
steps like Michael Shirts suggested.

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.

In the portion of the error below atom 9 is -C9-(H67,H68,H66) which in
this specific case H67 is already a dummy molecule with no mass or
charge.  From what I can see, it seems that the atoms do not know how
to treat the dummy molecules in terms of angles.  How should I treat
the dummy molecules? Should I be treating them like hollow spheres
with no charge so I would assign them angle constraints?

I think it can also be that I'm peturbing the dihedral angles
incorrectly.  I received errors at first saying that dihedral
multiplicities can't be peturbed so I had to equal the multiplicities
just to get it to run.  Does anybody have any experience with this?

Thank you for your help.

-Fabian Casteblanco

Portion of Error Output:
-----------------------------------------------------------------
Reading file nvt0.5.tpr, VERSION 4.5.3 (single precision)
starting mdrun 'SIMVASTATIN'
150000 steps,    300.0 ps.

Step 7, time 0.014 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.006111, max 0.139443 (between atoms 9 and 67)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length

Step 8, time 0.016 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.007341, max 0.167622 (between atoms 9 and 67)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length

Step 8, time 0.016 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.008771, max 0.201182 (between atoms 9 and 67)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length







On Mon, Oct 10, 2011 at 4:10 PM, Fabian Casteblanco
<fabian.casteblanco at gmail.com> wrote:
> Hi all,
>
> I have an additional question related to what Steven Neumann was
> mentioning.  I actually have to do a molecule mutation.  I'm trying to
> use Michael Shirts method  1) making small
> changes 'alchemical' changes in the molecules and computing the free
> energies by any method (BAR, TI, etc).  I'm specifically want to try
> to use BAR at the end once I collect all the data.  This helped a lot
> on clarification since it seemed that Justin's tutorial is essentially
> a FEP except its using the BAR mathematical method for computing the
> complete decoupling of the molecule rather than using the old FEP
> mathematics of the exponential averaging formula.  So BAR is only
> referring to the mathematical code used to calculate the overall free
> energy for the FEP, correct?
>
> 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)
>
> 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
>
>
> **Please, if anybody can help, I would greatly appreciate it.  Thanks.
> --
> Best regards,
>
> Fabian F. Casteblanco
> Rutgers University --
> Chemical Engineering
> C: +908 917 0723
> E:  fabian.casteblanco at gmail.com
>



-- 
Best regards,

Fabian F. Casteblanco
Rutgers University --
Chemical Engineering PhD Student
C: +908 917 0723
E:  fabian.casteblanco at gmail.com