AMBER Archive (2005)

Subject: Re: AMBER: A problem in TI approach - test case: meth2eth and eth2meth

From: Ilyas Yildirim (yildirim_at_pas.rochester.edu)
Date: Thu Apr 28 2005 - 19:00:11 CDT

Dear Prof. Case,

My problem is that the results depends on the cutoffs I chose in the
sander input file. Here is what I mean by that:

For both cases, I have used the following input file:

-------------------------------------------------------
    0 ps simulation with lambda = 0 - Test Run
     &cntrl
      imin=0,
      ntx=1,irest=0,
      ntpr=0,ntwr=0,ntwx=0,
      ntc=1,ntf=1,ntb=0,cut=30,
      igb=0,
      ntr=0,
      nstlim=0,dt=0.001,nscm=5000,nrespa=1,
      ntt=3,gamma_ln=1,tempi=300,temp0=300,ig=233,
      ntp=0,taup=1,pres0=1,
      icfe=1,clambda=0,klambda=1
     /
-------------------------------------------------------
In this particular case, the cut=30, and the system is in a non-periodic
boundary condition. And the results for eth->meth and meth->eth are what
I expected. Here are the results:

eth->meth -------------------------------------------------------
NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 331.57 PRESS = 0.0
Etot = 23.3991 EKtot = 26.6849 EPtot = -3.2858
BOND = 0.0029 ANGLE = 0.0007 DIHED = 0.0048
1-4 NB = -0.0210 1-4 EEL = 0.0109 VDWAALS = -1.5252
EELEC = -1.7590 EHBOND = 0.0000 RESTRAINT = 0.0000
DV/DL = 65.2824
------------------------------------------------------------------------------
meth->eth -------------------------------------------------------
NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 331.57 PRESS = 0.0
Etot = 23.3991 EKtot = 26.6849 EPtot = -3.2858
BOND = 0.0029 ANGLE = 0.0007 DIHED = 0.0048
1-4 NB = -0.0210 1-4 EEL = 0.0109 VDWAALS = -1.5252
EELEC = -1.7590 EHBOND = 0.0000 RESTRAINT = 0.0000
DV/DL = -65.2824
------------------------------------------------------------------------------

But the problems arises when I set the cut=10. Here are the results for
the case when cut=10:

eth->meth -------------------------------------------------------
NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 331.57 PRESS = 0.0
Etot = 23.3991 EKtot = 26.6849 EPtot = -3.2858
BOND = 0.0029 ANGLE = 0.0007 DIHED = 0.0048
1-4 NB = -0.0210 1-4 EEL = 0.0109 VDWAALS = -1.5252
EELEC = -1.7590 EHBOND = 0.0000 RESTRAINT = 0.0000
DV/DL = 65.2824
------------------------------------------------------------------------------
meth->eth -------------------------------------------------------
NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 331.57 PRESS = 0.0
Etot = 23.6633 EKtot = 26.6849 EPtot = -3.0215
BOND = 0.0029 ANGLE = 0.0007 DIHED = 0.0048
1-4 NB = -0.0210 1-4 EEL = 0.0109 VDWAALS = -1.2609
EELEC = -1.7590 EHBOND = 0.0000 RESTRAINT = 0.0000
DV/DL = -65.0181
------------------------------------------------------------------------------
As it can be seen, in meth->eth case, when cut=10, the results are not
same as the one found when cut=30. I was wondering what might be creating
this problem. Any idea is greatly appreciated. As far as I see from these
output results, the difference is coming from Van der Waals interactions,
but I dont know why this is happening. How can cut=10 play some role on
the differences?

All these simulations were done with the original 'sander'. I will do some
simulations with the 'modified sander' on some small structures, and will
use dummy atoms to annihilate/create hydrogen atoms. I will inform you all
with the results soon.

Best,

On Thu, 28 Apr 2005, David Case wrote:

> On Thu, Apr 28, 2005, Ilyas Yildirim wrote:
> >
> > The previous problem I have described was not about the modification I
> > have done by changing the source code. I just wanted to see the
> > initial (nstep=0) energy and dv/dl results of eth->meth and meth->eth
> > transformations. Theoretically, the dv/dl results for both cases should be
> > equal to each other for the same lambda value (with a sign difference),
> > but I have showed that this depends on the cutoff we chose. (In my case,
> > I have chosen 2 cutoffs, cut=10 and cut=30. In cut=30, the expected
> > result was satisfied, but for cut=10, it was not; which I could not
> > explain why) But in all these test cases, I used the original sander.
>
> You might be able to debug this in this way: Run short calculations at both
> lambda=1 and lambda=0 with both prmtop files. The lambda=1 results from one
> prmtop should exactly match the lambda=0 results from the other. If there are
> any differences, you will have to track them down to individual interactions
> that are present in one state but not in the other. You might have to print
> out all nonbonded pairs involving the solute atoms; fortunately, this is a
> pretty small list, so it seems likely that you could find out where the
> detailed differences are. Unlike Thomas, you seem(?) to be using a
> non-periodic calculation, so there may be no correspondence in the source of
> the problem.
>
> ...good luck...dac
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu
>
> 

-- 
  Ilyas Yildirim
  ---------------------------------------------------------------
  - Department of Chemisty       -				-
  - University of Rochester      -				-
  - Hutchison Hall, # B10        -				-
  - Rochester, NY 14627-0216     - Ph.:(585) 275 67 66 (Office)	-
  - http://www.pas.rochester.edu/~yildirim/			-
  ---------------------------------------------------------------

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu