Dear Ross,
thank you for the fast and exhaustive (as much as possible, due the
prototype nature of the NEB implementation) answer.
>
>> 1) after a careful reading of both AMBER9 manuals and Ross Walker's
>> Tutorial 7 from "Amber Workshop Tutorials" I understood that
>> NEB approach
>> can be performed using multiple frames from an approximate path as
>> starting guess for min. energy paths. In the case of our simulated
>> molecule, we obtained the desired conformational transition during a
>> standard MD simulation using implicit GB/SA solvation. Now,
>> we'd like to
>> obtain a NEB path using frames from this MD simulation. What
>> MD protocol
>> (in particular, temperature vs. time profile and frequency of frame
>> sampling from original MD run) is suggested in this case?
>
> This is the million dollar question. The problem is that the simulated
> annealing version of NEB in Amber 9 is new and because of this should be
> considered still a prototype method. That is I don't think we have a full
> understanding yet of the sampling issues involved or what is the best
> protocol to use. I don't think you can do much damage by over sampling in
> your simulated annealing but under sampling (running the simulated
> annealing
> for too short a time) will cause problems. Where the balance between this
> lies is very much dependent on your system size, degrees of freedom,
> number
> of images etc etc. Suffice to say that the larger the system is and the
> more
> images you have the longer you will need to run for. How long though would
> only be a guess.
>
> I would suggest you start your system at 0K and then do something along
> the
> follow lines. Note the large values of gamma_ln (this is required for NEB)
> to give the system viscosity. For a regular MD simulation these values are
> orders of magnitude too large:
>
> Step 1 -> Heat linearly from 0 to 300K dt=0.0005, skmin=10, skmax=10,
> gamma_ln=100,
> nstlim = 40,000 (20ps)
> Step 2 -> Heat linearly to 1000K, dt=0.002, skmin=20, skmax=20,
> gamma_ln=1000,
> nstlim = 25,000 (50ps)
> Step 3 -> Equilibrate at 1000K, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 4 -> Cool 1000K to 300K linearly, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 5 -> Equilibrate at 300K, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 6 -> Cool 300K -> 0K linearly, dt=0.002, skmin=25, skmax=25,
> gamma_ln=1000,
> nstlim = 50,000 (100ps)
> Step 7 -> 0K, dt=0.002, skmin=25, skmax=25, gamma_ln=5000, nstlim = 50,000
> (100ps),
> vv=1, vfac=0.1, ntwx=25000
>
> I have had success with this protocol. Although the time lengths are
> really
> just a guess. You may want to increase the length of sampling in the
> middle
> if you can afford it. You should also probably run 10 or so independent
> runs
> with different random seeds to test how well your sampling is converged.
> You
> can then see how many different reactions paths you get and compare them.
>
I'll try your protocol, possibly with systematic exploration of some of
the involved parameters, and, if you are interested, I can send you the
main results/issues found in the analysis.
> In terms of building your initial guess at the pathway there ar eno hard
> and
> fast rules. You just have to provide the code with sufficient input
> structures. And indeed this is very much a "guess" at the pathway so the
> structures really just need to be reasonable representations of points
> along
> the pathway. So I would just pick what you thing are representative of the
> transition from your MD simulation and string them together.
I just wondered if some "quasi-minimization-only" protocol could be
meaningful for opportunely-sampled structures from a plain MD simulation
involving the transition under examination.
> You may also want to just try it with half the images piled up on one end
> and half on the other end (I.e. no guess at the pathway). We have had some
> success with this. You obviously need to equilibrate for longer but if you
> are lucky it won't get trapped in a local minimum due to the high
> temperature.
>
Calculations with this approach are already running, as we were interested
in a comparison between the two approaches and for some of the systems
under examination we haven't the MD trajectory including the transition.
>> 2) looking at the source code of AMBER9 implementation of the quenched
>> velocity Verlet algorithm used for NEB minimization (see below), we
>> observed that, differently from what stated in the manual,
>> when the dot
>> product v dot f<0 , velocities are simply zeroed and not
>> scaled according
>> to eq. (6.34) v=x(v dot f)f and, consequently, the input
>> scaling factor
>> VFACT is never used (or, equivalently, the program behaves as if VFACT
>> value were hardwired to 0). I've found no explanation for
>> this discrepance
>> either in docs, or in man. upgrades, or in mailing list, or
>> in bugfixes.
>>
>> >>>> File: runmd.f
>> >>>> Subroutine: quench
>> > if (dotproduct>0.0d0) then
>> > v(1:3*natom) = dotproduct*f(1:3*natom)*force
>> > else
>> > v(1:3*natom) = 0.0d0
>> > end if
>
> This is a typo in the manual. It should state:
>
> if(v.f>=0) then v = (v.f)f
>
> if(v.f<0) then v = 0
>
> This is what the code does and is what should happen. I will update the
> manual errata page.
OK
>
> Note you may also want to read the following paper describing the NEB
> implementation in Amber 9:
>
> Mathews, D.H., Case, D.A. J. Mol. Biol. (2006) 357, 1683-1693
>
We already read the paper and most of the information was desumed from it.
Thank you again,
sincerely,
Pietro
--
Dr. Pietro Amodeo
Istituto di Chimica Biomolecolare del CNR
Comprensorio "A. Olivetti", Edificio 70
Via Campi Flegrei 34
I-80078 Pozzuoli (Napoli) - Italy
Phone +39-0818675072
Fax +39-0818041770
Email pamodeo.icmib.na.cnr.it
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Oct 01 2006 - 06:07:16 PDT
