[gmx-users] Peturbing a Dihedral for FEP

[Fabian Casteblanco]

Hello all,

 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
>



-- 
Best regards,

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