Re: [AMBER] TI simulations blowing up with clambda = 1 with pmemd

From: Stefan Ivanov <stefan.ivanov.postgrad.manchester.ac.uk>
Date: Sun, 29 Jan 2017 15:07:13 +0000

Hi Hannes, hello all,

The trajectories seem to be running fine now, but at lambda = 1, during the first few thousand steps of heating, the MBAR energies for 0.0 skyrocket:

step 1000:
MBAR Energy analysis:
Energy at 0.0000 = -28599.9657
Energy at 0.1000 = -28795.7667
Energy at 0.2000 = -28843.1242
Energy at 0.3000 = -28862.9072
Energy at 0.4000 = -28873.4368
Energy at 0.5000 = -28879.8249
Energy at 0.6000 = -28884.0013
Energy at 0.7000 = -28886.8353
Energy at 0.8000 = -28888.7562
Energy at 0.9000 = -28889.9686
Energy at 1.0000 = -28890.5172


Step 2000:

MBAR Energy analysis:
Energy at 0.0000 = -27983.8250
Energy at 0.1000 = -28271.9302
Energy at 0.2000 = -28336.4348
Energy at 0.3000 = -28365.3598
Energy at 0.4000 = -28381.8330
Energy at 0.5000 = -28392.4822
Energy at 0.6000 = -28399.9612
Energy at 0.7000 = -28405.5528
Energy at 0.8000 = -28409.9725
Energy at 0.9000 = -28413.7016
Energy at 1.0000 = -28417.2183


Step 3000:

MBAR Energy analysis:
Energy at 0.0000 = -18780.0639
Energy at 0.1000 = -27628.3326
Energy at 0.2000 = -28036.1619
Energy at 0.3000 = -28157.6530
Energy at 0.4000 = -28211.9390
Energy at 0.5000 = -28241.3843
Energy at 0.6000 = -28259.3539
Energy at 0.7000 = -28271.2685
Energy at 0.8000 = -28279.7113
Energy at 0.9000 = -28286.0920
Energy at 1.0000 = -28291.3495

Step 4000:

MBAR Energy analysis:
Energy at 0.0000 = 1397837.9002
Energy at 0.1000 = -26377.2229
Energy at 0.2000 = -27601.2916
Energy at 0.3000 = -27901.9458
Energy at 0.4000 = -28021.5126
Energy at 0.5000 = -28080.7553
Energy at 0.6000 = -28114.2043
Energy at 0.7000 = -28134.8575
Energy at 0.8000 = -28148.5085
Energy at 0.9000 = -28158.0694
Energy at 1.0000 = -28165.1784

Step 5000:



MBAR Energy analysis:
Energy at 0.0000 = 274909.2086
Energy at 0.1000 = -26760.0263
Energy at 0.2000 = -27654.5419
Energy at 0.3000 = -27884.5697
Energy at 0.4000 = -27975.7088
Energy at 0.5000 = -28019.8213
Energy at 0.6000 = -28043.7215
Energy at 0.7000 = -28057.5592
Energy at 0.8000 = -28065.8453
Energy at 0.9000 = -28070.8135
Energy at 1.0000 = -28073.6622

Step 6000:

MBAR Energy analysis:
Energy at 0.0000 = ************
Energy at 0.1000 = -26265.1456
Energy at 0.2000 = -27509.2723
Energy at 0.3000 = -27811.5121
Energy at 0.4000 = -27930.5115
Energy at 0.5000 = -27988.1462
Energy at 0.6000 = -28019.2967
Energy at 0.7000 = -28037.1517
Energy at 0.8000 = -28047.6104
Energy at 0.9000 = -28053.6329
Energy at 1.0000 = -28056.8779

The extremely high values carry over into the subsequent density equilibration and production simulations. Conversely, for low lambda values, the high values are at 1.0:

Lambda = 0, step 6000:

MBAR Energy analysis:
Energy at 0.0000 = -27959.9269
Energy at 0.1000 = -27954.6351
Energy at 0.2000 = -27948.3082
Energy at 0.3000 = -27940.4340
Energy at 0.4000 = -27930.3231
Energy at 0.5000 = -27916.5964
Energy at 0.6000 = -27896.0779
Energy at 0.7000 = -27860.1191
Energy at 0.8000 = -27777.1385
Energy at 0.9000 = -27432.2432
Energy at 1.0000 = 3381772.1010

Step 17000:

MBAR Energy analysis:
Energy at 0.0000 = -27646.8667
Energy at 0.1000 = -27642.2782
Energy at 0.2000 = -27636.0701
Energy at 0.3000 = -27627.7974
Energy at 0.4000 = -27616.5953
Energy at 0.5000 = -27600.5964
Energy at 0.6000 = -27575.4905
Energy at 0.7000 = -27529.7575
Energy at 0.8000 = -27422.8804
Energy at 0.9000 = -26997.0203
Energy at 1.0000 = ************

Although the trajectories seem fine, this certainly doesn't look like sensible output. I tried running with and without SHAKE, in NVT, and in NPT, but the picture is pretty much the same.

My input files:

minimization:

&cntrl
imin = 1, ntmin = 2, maxcyc = 1000, ncyc = 200,
ntwe = 100, ntpr = 10, ntwr = 100,
ntc = 2, ntf = 1,
noshakemask = ':1XJ,1X3',
icfe = 1, clambda = 1.0, scalpha = 0.5, scbeta = 12.0,
timask1 = ':1XJ', timask2 = ':1X3',
ntr = 1, restraint_wt = 999.00,
restraintmask=':1XJ,1X3',
ifsc=1,
scmask1=':1XJ.C4,C5,C14,C15,C17,C24,C25,C26,C28,C29,C30,C37,C38,C39,C40,C41,C42,C47,N49,O58,S63,Cl6,F59,F60,F61,H4,H5,H14,H15,H49,H281,H282,H291,H292,H301,H302,H371,H372,H381,H382,H401,H402,H403,H411,H412,H413,H421,H422',
scmask2=':1X3.N1,N93,O1,C81,C82,C83,C4,C5,C86,C87,C88,C89,C90,H35,H36,H81,H82,H83,H4,H5,H86,H87,H88,H89',
crgmask='',
/


heating:

heating
 &cntrl
   imin = 0, nstlim = 100000, irest = 0, ntx = 1, dt = 0.001,
   ntt = 1, temp0 = 300.0, tempi = 50.0, tautp = 1.0,
   ntc = 2, ntf = 1, noshakemask=':1XJ,1X3', ig = -1,
   ntb = 1,
   ioutfm = 1, iwrap = 1,
   ntwe = 1000, ntwx = 1000, ntpr = 1000, ntwr = 5000,
   icfe = 1, clambda = 1.0, scalpha = 0.5, scbeta = 12.0,
   ifmbar = 1, bar_intervall = 1000, bar_l_min = 0.0, bar_l_max = 1.0,
   logdvdl = 1,
   timask1 = ':1XJ', timask2 = ':1X3',
   ntr = 1, restraint_wt = 999.00,
   restraintmask=':1XJ,1X3',
   ifsc=1,
   scmask1=':1XJ.C4,C5,C14,C15,C17,C24,C25,C26,C28,C29,C30,C37,C38,C39,C40,C41,C42,C47,N49,O58,S63,Cl6,F59,F60,F61,H4,H5,H14,H15,H49,H281,H282,H291,H292,H301,H302,H371,H372,H381,H382,H401,H402,H403,H411,H412,H413,H421,H422',
   scmask2=':1X3.N1,N93,O1,C81,C82,C83,C4,C5,C86,C87,C88,C89,C90,H35,H36,H81,H82,H83,H4,H5,H86,H87,H88,H89',
   crgmask='',
 /

 &ewald
 /

 &wt
   type='TEMP0',
   istep1 = 0, istep2 = 100000,
   value1 = 50.0, value2 = 300.0
 /

 &wt type = 'END'
 /


Any advice would be immensely appreciated.

Sincerely,

Stefan


________________________________________
From: Hannes Loeffler [Hannes.Loeffler.stfc.ac.uk]
Sent: Friday, January 27, 2017 3:46 PM
To: amber.ambermd.org
Subject: Re: [AMBER] TI simulations blowing up with clambda = 1 with pmemd

On Fri, 27 Jan 2017 14:54:38 +0000
Stefan Ivanov <stefan.ivanov.postgrad.manchester.ac.uk> wrote:

> Hi Hannes,
>
> The previous e-mail was too big to post on the mailing list; it
> wouldn't let me.

Ok.


> Hmmm, the minimization seems to be doing the trick, although there
> aren't any clashes in my starting structures. Also, at lamda = 1 BOTH
> ligands would become completely distorted although at that lambda
> value the template isn't interacting with its environment at all.
> Also, even if there were clashes between the modeled ligand and any
> water molecules and/or ions, it wouldn't cause the whole molecule to
> become completely distorted and keep twisting and deforming
> throughout the whole run, let alone affect the template ligand at
> lambda = 1.

Large parts of your two ligands are single coordinate (or effective
single topology if you like) and you have both appearing and
disappearing parts (softcore). So most of your molecules are present
at all times and single coordinates will distort together.

I have only checked for 20000 steps of heating and 50000 steps of
pressurising. I don't see anything that I would call distorted. One
thought: what if instead of the lib files you have created you directly
load the mol2 files instead and/or update to AmberTools16? BTW, I am
on Amber16.


> I've run 2000 steps of minimization on similar structures and I've
> had similar problems so I'm not sure if this is a general solution or
> works just in this particular case. I'm not really sure what's going
> on here.

I think all you need is a few tens of minimisation. When I set up the
system I had several waters <= 2A near the ligands.

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sun Jan 29 2017 - 07:30:02 PST
Custom Search