Bizarre behavior in QMMM optimization


Clicked A Few Times
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

Clicked A Few Times
Hard to say what went wrong. Why don't you send me all the relevant files so that I can reproduce it on my end.

Marat

Clicked A Few Times
Sure thing. I just e-mailed them to you. Thanks!!

Clicked A Few Times
The problem was in your amber.par file. The line referring to HX was shifted by one space to the right. The classical MD module unfortunately assumes a rigid format for all the files involved and it is VERY IMPORTANT to adhere to that.

Also I would recommend to specify cutoff distance in the MD block, which should less than half of your box size, e.g.


md
system amine_opt
cutoff 0.5 qmmm 0.5
end

Marat

Clicked A Few Times
Thanks Marat! I'm really not sure how I missed that...


Forum >> NWChem's corner >> QMMM