Dear Amber Mailing List,
I have done some perturbation simulations, and realized something strange
in the results. I will first describe the system, how I created them, and
then describe the problem.
I have done 2 different perturbations; one from ethane --> methane and the
other one from methane --> ethane. I know that 'sander' isnt gonna work
when I have dummy atoms in the initial state, but my md simulations are
'zero' step simulations. I am doing these simulations to test some
theoretical expectation; like if I am gonna get the same free energy
difference with different signs for both cases, or if the total energy of
one perturbation at lambda=0.5 is equal to the total energy of the other
perturbation at lambda=0.5, etc.
First, I have created an eth2meth.pdb file using xleap. eth2meth.pdb file
has an ethane molecule and some water molecules around it. I am attaching
this file to my email. (eth2meth.pdb)
Then, I opened this file in xleap, with frcmod.dummy file. (frcmod.dummy
has the dummy atom parameters in it. I am attaching this file with another
file, eth.lib, which is the library for my eth residue.) Then I created
the eth2meth.prmtop and eth2meth.inpcrd files. I have created the
meth2eth.prmtop and meth2eth.inpcrd files starting from this initial
eth2meth.pdb file. I have manipulated the structure such that initially it
has the methane molecule and will be transformed to an ethane molecule. (I
am attaching the leap.log, the eth2meth.prmtop, eth2meth.inpcrd,
meth2eth.prmtop and meth2eth.inpcrd files)
Finally, I have created exactly the same .inpcrd files for both eth2meth
and meth2eth simulations. The reason why I wanted to have exactly the same
inpcrd for both perturbations is to be compatible when comparing the
results. (Here, someone has to be careful. In eth2meth perturbation, the
initial structures does not have any dummy atoms in it, but at the final
state it has it. The -C-H3 will be transformed into -H-(DH)3 where DH
stands for the dummy atom. And for the meth2eth perturbation, the initial
structure has 3 dummy atoms in it, but at the final state, they will be
transformed to -C-H3. But in both perturbations, the initial coordinates
are the same. In case it is not clear, use the prmtop and inpcrd files to
create the pdb structure for the molecules, such as,
ambpdb -p meth2eth.prmtop <meth2eth.inpcrd> meth2eth.pdb
and
ambpdb -p eth2meth.prmtop <eth2meth.inpcrd> eth2meth.pdb
)
Then, for both perturbations, eth2meth and meth2eth, I have used 64
clambda values. I am attaching the script I used to create these
directories/files. (create_files.pl)
For all of the simulations, I have used the following test.in input file
for the sander simulation:
------------------ test.in----------------
0 ps simulation with lambda = $lambda - Test Run
&cntrl
imin=0,
ntx=1,irest=0,
ntpr=0,ntwr=0,ntwx=0,
ntc=1,ntf=1,ntb=0,cut=10,
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=$lambda,klambda=6
/
------------------------------------------
Here, $clambda is the lambda values I have chosen. I am including the
script file (test_job.pl) which creates these input files for each clambda
value, and runs the sander job.
As it can be seen, this is a zero step simulation (nstlim=0). Why did I do
a simulation like this? Theoretically, I would expect the following
results tobe seen:
1. free energy diff(meth-->eth) = - free energy diff(eth-->meth) [After
using trapezoidal rule]
2. (for each lambda value) dvdl(meth-->eth) = - dvdl(eth-->meth)
3. (for k=1 and lambda=0.5) E_total(meth-->eth) = E_total(eth-->meth)
Here are my results:
------------------------------------------
(When 'cut' in test.in was set to 10)
1. Didnt satisfied (see meth2eth_lambda_vs_dvdl_cut_10 and
eth2meth_lambda_vs_dvdl_cut_10 files)
2. Didnt satisfied (see meth2eth_lambda_vs_dvdl_cut_10 and
eth2meth_lambda_vs_dvdl_cut_10 files)
3. Didnt satisfied. (I did not include the files here, but it is a very
fast simulation)
------------------------------------------
(When 'cut' in test.in was set to 30)
1. Was satisfied (see meth2eth_lambda_vs_dvdl_cut_30 and
eth2meth_lambda_vs_dvdl_cut_30 files)
2. Was satisfied (see meth2eth_lambda_vs_dvdl_cut_30 and
eth2meth_lambda_vs_dvdl_cut_30 files)
3. Was satisfied (Again, I did not include the files here, but it is a
very fast simulation)
------------------------------------------
So, the question is the following:
Given that we have 2 perturbatios (eth-->meth and meth-->eth) and that we
do not use any boundaries and that we have the same .inpcrd files for both
perturbations, why arent the 3 expectations satisfied when cut=10? I could
not get any answer for these questions, yet. It is interesting that when
cut=30, then I dont have any problem.
The main problem is when I use a periodic boundary in my system. There is
a limit on the 'cut' in a periodic boundary simulation. But if the results
depends on the 'cut' (as it happens in the simulations I have done), I am
not sure how to overcome this problem.
Thanks in advance for your answers. Best,
(PS: I could not send eth.lib file, because it was always giving an
errors, saying that
<amber.scripps.edu>: host relay1.scripps.edu[137.131.200.29] said:
554 5.7.1
E-mail attachment had bad filename 'eth.lib' (in reply to end of DATA
command)
)
--
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.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Tue Apr 26 2005 - 02:53:00 PDT
