AMBER Archive (2007)

Subject: Re: AMBER: SCF convergence issues

From: Steven Winfield (saw44_at_cam.ac.uk)
Date: Fri Feb 23 2007 - 13:13:33 CST

Hi Ross,

Thanks for your reply.

> A cutoff of 6.5 angstroms for MM is dangerously small and a cutoff of 3.5
> angstroms for the qm is very very small. I would never reduce this below 8
> angstroms. I know this can give problems with distributed QM regions in
> small periodic boxes but it is there for a reason. The cutoff cannot extend
> beyond half the box or you get periodicity artifcacts and you can't make it
> less than 8 angstroms or you seriously truncate the interactions. Thus the
> only really acceptable solution is to make the box bigger.

I ran some tests on a purely classical system, with the cutoff
increasing from 3.0A to around 12.0A and found that by about 6.0A there
were very little changes in the forces. My understanding of this
variable was that it just divides the coulombic interactions between
real and reciprocal space. If so, then shouldn't changing this variable
just alter the speed at which the PME is calculated?

> Most of the scf cycles converge properly. It is only one of them that don't.
> However, the difference in energy between two steps is huge. E.g.
> QMMM: SCF Energy = -1325.32549345000734 KCal/mol, -5545.16186459483106
> KJ/mol
> QMMM: SCF Energy = -1298.86881677272322 KCal/mol, -5434.46712937707434
> KJ/mol

Maybe you've answered this in other way, but why are there multiple sets
of SCF iterations? The input files specify there should be only one
integration step, so for there to be multiple SCF iterations AFTER the
debug force routine seems to have dumped forces seems confusing.

> Is this system minimized? It looks to me that this system may be a long way
> from equilibrium and that is causing the problems. Did you minimize it (even
> 100 steps of steepest descent would be good). I would then try some simple
> MD and see how it behaves.

I built the system from a few TIP3PBOXs, manually adjusted the box
dimensions to give the correct density at 300K, performed a few thousand
minimisation steps (about half steepest descent), then performed MD from
0K -> 300K over 100ps, followed by another 100ps at 300K, both NVT.

> The above problem could also be caused by atoms moving in and out of the
> cutoff. With such a small cutoff you will get huge jumps in energy as an
> atom crosses the cutoff and this will lead to instabilities.

But this is one integration step... no atoms are moving.

> I would also try turning pseudo_diag back on as this can often get you away
> from bad initial densities quickly. The issue above is that for the SCF set
> that doesn't converge the initial density matrix is very bad. This is likely
> caused by the system moving a long way from the previous step. I.e. several
> atoms moved outside the small cutoff causing a large change in the field and
> the difference was so severe that the scf routine couldn't recover.

I experenced the same problems with pseudo_diag on, and the default
scfconv, and itrmax=10000 (and combinations of these)

> On another note I don't think I ever tested dumpfrc with QM/MM. This tries
> to decompose the energies and forces and this is not possible with QM/MM,
> since you can't pairwise decompose things and so this could also be the
> route of your problem. Why are you trying to decompose the force array? If
> you want to test the accuracy of the forces you can use do_debugf=1 and
> atomn=1,2,3 etc. You probably want to set scf_conv=1.0D-8 and tight_p_conv=1
> though to obtain decent gradients. 10-6 is probably not tight enough. I'd
> also just try some simple MD first.

I'm involved in the writing of program (or rather a library) which
implements a new QM/MM integration scheme. Part of this scheme is the
ability to use existing programs as 'black boxes' for force calculations.

Currently, my program writes an AMBER input file and calls sander with
the dumpfrc option, then parses the forcedump.dat files for the forces
on each atom. In that respect I don't mind if the decomposition is
correct, since I only care about the total force on each atom which
AMBER itself would use for integrating the equations of motion (and
which I also assumed the ones in the relevant section of forcedump.dat).

If there is a 'safer' or quicker way of obtaining the forces for a given
input configuration (without modifying the code) then I would very much
like to know.

Best regards,

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