AMBER Archive (2007)

Subject: RE: AMBER: SCF convergence issues

From: Ross Walker (ross_at_rosswalker.co.uk)
Date: Fri Feb 23 2007 - 17:35:47 CST

Dear Steven,

> 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?

Right, but this also controls the van der Waals cutoff since the VDW and
direct space EEL interactions are evaluated at the same time. Thus you are
truncating the VDW interactions at 6.0 angstroms which is not the best idea.
If you are just testing things in water you are likely seeing skewed results
here because TIP3P etc have no VDW radius on the hydrogens. They are
parameterized this way. Thus a 6 angstrom cutoff of VDW for pure water isn't
strictly a 6 angstrom cutoff since it is only the distance between oxygens
that matter here. However, if you later switch to a protein then you will
see issues. Hence 8 angstroms is normally the smallest recommended in terms
of energy conservation etc. At 6 angstroms for a protein in water I think
you will see energy conservation errors from the truncation of the VDW
forces.

> > 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?

Right, exactly, I don't know... The only portion of the debug code I have
ever tested with QM/MM is the numerical force comparisson. Here numerous SCF
calculations are done as each atom whose force is being tested is moved a
small distance in the x,y and z directions. I have never looked at the
dumpfrc routines so it is possible that this does some kind of seperate
calculation, e.g. runs the bond calculation, then the angles etc etc. I'd
have to have a good look at this code to work out what is going on but I
suspect as things are right now this section of code is probably not
compatible with the QM/MM code. I'll have to try it out myself when I get
some time. I'd like to know if you can run simple MD with your system though
as this will illustrate if it is a problem with the debug code and not just
something strange with your system that is causing covergence problems. If
everything else works okay then I'll try and see if there is an easy way of
getting dumpfrc to work.

>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.

Indeed, there is something not right here... It is obviously not just
dumping the frc array but re-calling parts of the code, including the QM/MM
routines several times.

> 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.

Was this with QM/MM MD? If so then is definately a problem with
incompatibility ofr dumpfrc and QMMM.

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

Yes I see now. I initially thought that it was also trying to do numerical
gradients by moving each atom in turn and comparing the energy.
 
> I experenced the same problems with pseudo_diag on, and the default
> scfconv, and itrmax=10000 (and combinations of these)

I think by the time you get to the scf that doesn't converge everything is
already corrupted at this point. I wouldn't be surprised if the initial
density matrix guess is just completely scrambled by the debug code calling
things in strange orders etc. This is especially true if you only see this
when you are in debug mode.
 
> 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.

The safer way would just be to put a write(20,*) f(1:3*natom) at the end of
the force routine followed by a call to mexit but this obviously involved
adding a couple of lines to force.f. I don't know if you would be okay to do
this or not. Failing that I can try and fix the dumpfrc option but it could
be quite a bit of work and would take a while before I have a chance to look
at it.

All the best
Ross

/\
\/
|\oss Walker

| HPC Consultant and Staff Scientist |
| San Diego Supercomputer Center |
| Tel: +1 858 822 0854 | EMail:- ross_at_rosswalker.co.uk |
| http://www.rosswalker.co.uk | PGP Key available on request |

Note: Electronic Mail is not secure, has no guarantee of delivery, may not
be read every day, and should not be used for urgent or sensitive issues.

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