Hi everyone...
I've run into some very strange behavior in QMMM optimization that I can reproduce in a smallish molecule; hopefully someone has seen similar behavior or can help me figure out what's going on.
I start with a PDB file for CH3NH2 (amine.pdb):
ATOM 1 CB1 UNK X 1 -3.336 -0.266 -2.632 1.00 0.00
ATOM 2 HB1 UNK X 1 -4.313 -0.467 -2.179 1.00 0.00
ATOM 3 HB2 UNK X 1 -3.403 -0.526 -3.703 1.00 0.00
ATOM 4 HB3 UNK X 1 -2.616 -0.950 -2.168 1.00 0.00
ATOM 5 NB1 UNK X 1 -2.950 1.119 -2.344 1.00 0.00
ATOM 6 HB4 UNK X 1 -2.049 1.319 -2.777 1.00 0.00
ATOM 7 HB5 UNK X 1 -3.619 1.767 -2.763 1.00 0.00
END
and prepare it with the following script:
start amine
prepare
source amine.pdb
new_top new_seq
new_rst
modify atom 1:_CB1 quantum
modify atom 1:_HB1 quantum
modify atom 1:_HB2 quantum
modify atom 1:_HB3 quantum
modify atom 1:_NB1 quantum
modify atom 1:_HB4 quantum
modify atom 1:_HB5 quantum
center
orient
solvate box 1.0
update lists
ignore
write amine_ref.pdb
write amine_ref.rst
end
task prepare
I followed the procedure described in Marat's tutorial, creating UNK.frg:
# This is an automatically generated fragment file
# Atom types and connectivity were derived from coordinates
# Atomic partial charges are crude estimates
# 12/05/11 17:37:17
#
$UNK
7 1 1 0
UNK
1 CB1 CX 0 0 0 1 1 -0.150000 0.000000
2 HB1 HX 0 0 0 1 1 0.050000 0.000000
3 HB2 HX 0 0 0 1 1 0.050000 0.000000
4 HB3 HX 0 0 0 1 1 0.050000 0.000000
5 NB1 NX 0 0 0 1 1 0.000000 0.000000
6 HB4 HX 0 0 0 1 1 0.000000 0.000000
7 HB5 HX 0 0 0 1 1 0.000000 0.000000
and amber.par:
Electrostatic 1-4 scaling factor 0.833333
Relative dielectric constant 1.000000
Parameters epsilon R*
#
Atoms
CX 12.01000 3.59824E-01 1.90800E-01 1 1111111111
Q 6 1.79912E-01 1.90800E-01
NX 14.01000 7.11280E-01 1.82400E-01 1 1111111111
Q 7 3.55640E-01 1.82400E-01
HX 1.00800 6.56888E-02 6.00000E-02 1 1111111111
Q 1 3.28444E-02 6.00000E-02
End
When I run my preparation script, all goes according to plan and amine_ref.pdb looks sensible. Now comes the interesting part. I want to optimize my solvent/solute system, so I run the following script (after copying amine_ref.rst to amine_opt.rst):
start amine_opt
permanent_dir ./perm
scratch_dir ./scratch
charge 0
md
system amine_opt
end
basis
C library 6-31G*
N library 6-31G*
H library 6-31G*
end
dft
xc b3lyp
direct
iterations 500
end
qmmm
region qm solvent
eatoms -94.3220090769
maxiter 50 5000
density espfit
ncycles 100
xyz amine_opt
end
task qmmm dft optimize
Everything looks peachy until I visualize the minimization trajectory saved in my XYZ files. Then I see that the molecule is minimizing into a really strange configuration, with hydrogens that have been bonded to the same atom collapsing to the same position. The last XYZ file I got before I gave up and killed the optimization was:
7
geometry
C -0.56864875 -0.80610217 -0.29876376
H 0.77935139 -0.73503754 1.18450010
H 0.77865950 -0.73247840 1.18417894
H 0.77848966 -0.73270980 1.18405198
N -1.14914162 0.02864728 0.29775924
H -0.50502373 1.22595179 0.34530901
H -0.50470915 1.22577115 0.34500492
I suspect that the problem is in the energy calculation, because when I look at the minimization progress using 'grep @' I see the following:
@------------------------------------------------
@
@ QM/MM Multiple Region Optimization
@ region1: qm with bfgs maxiter = 50
@ region2: solvent with sd maxiter =***
@
@
@------------------------------------------------
@ ncycle = 1
@
@
@ Optimizing qm region with bfgs
@
@ Step Energy Delta E Gmax Grms Xrms Xmax Walltime
@ ---- ---------------- -------- -------- -------- -------- -------- --------
@ 0***************** 0.0D+00****************** 0.00000 0.00000 4.3
@ 1***************** -1.1D+21****************** 0.31224 0.36907 20.0
@ 2***************** -2.6D+20****************** 0.30228 0.43258 23.4
@ 3***************** -2.1D+20****************** 0.26719 0.64393 31.6
@ 4***************** -9.3D+19****************** 0.24605 0.41148 37.1
... etc.
I also see those asterisks instead of energy output whenever other systems show this same problem, but following the same process (as far as I can tell) gives me normal energy output and sane-looking optimization for some systems (a CH4 molecule in a water box, for example).
Am I missing something really important here? Any tips at all would help a lot.
Thanks (as usual),
--craig
|