MM atom collapsing into QM nucleus in Tinker/Molcas

ArnoldBax · 2018-05-23 12:16:56

I'm running a QM/MM optimisation, using HF and Amber99 force field with Molcas7.8/Tinker5.1.
I am moving both the QM system and some MM atoms. All seems to be fine, but after some optimisation steps a hydrogen atom of a MM water molecule close to the QM system ends up on top of a QM carbon nuclei, and the calculation stops.

In Amber, polar hydrogens have zero VdW parameters so I guess the problem is that the electrostatic interaction between the MM water hydrogen and the QM carbon overcomes the penalty due to stretching of the O-H bond, also because the ESPF charge on the QM carbon gets large.

A simple solution would be to set non-zero VdW parameters on polar hydrogen atoms, though it may require some parameterization.
Has anyone had a similar problem, and found a way to avoid this?

Thanks!

Here's Molcas input. In the Tinker key file QMMM microiteration is ON.

  &Gateway 
Tinker 
Basis  =  6-31G*  
Group  =  Nosym 

>>>  Export  MOLCAS_MAXITER=200
>>>  Do  while 

  &SEWARD 
  &Espf ; External  =  Tinker ; LAMorok ; Show
  &SCF &END
  Title
  WQ 
  Charge
  -1 
  End of input
  &SLAPAF ; Cartesian

>>>  EndDo

niko · 2018-05-23 13:27:41

Hi,
I never met this kind of problem before. Do you include this particular water molecule into the list of MM atoms known by Molcas (ie using the QMMM keyword in Tinker's key file)? In that case, Molcas optimizer is taking care of it and may be faulty. If not, the problem is occurring during Tinker microiterations? Could you share your key file?

ArnoldBax · 2018-05-23 14:21:49

Thank you for the reply. I think the problem is occurring during Tinker microiterations. Here's the .key file.... I also tried with Steepest Descent but the problem is still there.

parameters /software/tinker/tinker-5.1.09-qmmm/params/amber99sb.prm
QMMM 30
QM -923 944 916 917 
MM  914 918 
MM  945
LA 15828 15829 15830
atom   2999    99    HLA     "Hydrogen Link Atom"        1      1.008     0
active -914 946 *move whole chro residue
active -15828 15830
active -146 199 
[... several other active indexes...]
active 13548 13549 13550
charge -914 0.0
charge -915 0.0
charge -918 0.0
charge -945 0.0
charge -946 0.0
charge -919 -0.0441
charge -920 -0.0441
charge -921 -0.5128
charge -922  0.4246
QMMM-MICROITERATION ON
QMMM-ELECTROSTATICS ESPF

VERBOSE

this is what Tinker prints at the last step (you can see the energies getting huge)

  Limited Memory BFGS Quasi-Newton Optimization :
 
  QN Iter    F Value      G RMS     F Move    X Move    Angle  FG Call  Comment
 
      0   -4518.2342     1.2356                                    1
      1   -4519.2650     1.2215     1.0308    0.0014     0.00      3    Success
      2   -4520.0957     0.9400     0.8306    0.0010    44.73      4    Success
      3   -4521.4638     1.0396     1.3681    0.0021    44.43      5    Success
      4   -4523.7999     1.7288     2.3362    0.0040    56.79      6    Success
      5   -5274.2449   614.2984   750.4449    0.0454    72.49     18    ReSearch
      6   -5274.2475   614.7376     0.0026    0.0000    89.95     20    Success
      7   -5374.3691   576.8959   100.1216    0.0014    83.79     28    Success
      8   -6565.5444  1001.6759  1191.1753    0.0072    71.62     32    Success
      9   -7714.1637  1410.2126  1148.6193    0.0084    76.28     33    Success
     10  -15459.6738  7963.9516  7745.5101    0.0164    78.25     36    Success
     11  -15766.9322  8406.8338   307.2584    0.0016    88.03     38    Success
     12  -32620.4159 43933.2668 16853.4837    0.0161    86.30     44    Success
     13  -33110.3089 45506.0143   489.8930    0.0013    89.22     46    Success
     14  -37091.1712 57330.2975  3980.8624    0.0041    88.75     48    Success
     15  -0.4995D+05 0.1054D+06 0.1286D+05    0.0189    89.03     50    Success
     16  -0.7787D+05 0.2659D+06 0.2792D+05    0.0170    89.24     53    Success
     17  -0.8397D+05 0.3117D+06 0.6096D+04    0.0045    89.59     55    Success
     18  -0.1409D+06 0.8874D+06 0.5695D+05    0.0054    88.44     59    Success
     19  -0.1440D+06 0.9266D+06 0.3094D+04    0.0011    89.21     61    Success
     20  -0.1485D+06 0.9848D+06 0.4464D+04    0.0008    89.61     62    Success
     21  -0.2956D+06 0.3918D+07 0.1471D+06    0.0216    89.74     67    Success
     22  -0.3051D+06 0.4173D+07 0.9461D+04    0.0031    89.94     69    Success
     23  -0.7388D+06 0.2477D+08 0.4338D+06    0.0250    89.88     73    Success
     24  -0.7490D+06 0.2546D+08 0.1020D+05    0.0016    89.98     75    Success
     25  -0.1312D+07 0.7831D+08 0.5630D+06    0.0020    89.52     81   SmallFct

ArnoldBax · 2018-05-23 18:06:53

I tried with a minimal system, 1 QM water and 1 MM water. With numerical gradients it works fine, while with analytical gradients (invoked if the Cholesky keyword is omitted) I get the same problem as before.

file.inp is

  &Gateway 
Tinker 
Basis  =  6-31G*  
Group  =  Nosym 

&SEWARD &END
 Title = Wat
 CHOlesky
End of input

  &Espf 
External  =  Tinker 
LAMorok 

>>>  Export  MOLCAS_MAXITER=200
>>>  Do  while 

  &SEWARD ; CHOlesky
  &Espf ; External  =  Tinker ; Show
  &SCF &END
  Title
  WQ 
  Charge
  0  
  End of input
  &SLAPAF ; Cartesian

>>>  EndDo

file.xyz

 6   molden generated tinker .xyz (amber param.)
 1  O     9.170000   9.040000   5.030000  2001   2  3
 2  H     9.500000   9.030000   4.120000  2002   1
 3  H     8.270000   8.710000   4.960000  2002   1
 4  O     6.490000   8.070000   5.130000  2001   5  6
 5  H     6.440000   7.270000   5.670000  2002   4
 6  H     5.790000   8.620000   5.470000  2002   4

file .key

QMMM 3 
QM -4 6 
active -1 6 
QMMM-MICROITERATION ON
QMMM-ELECTROSTATICS ESPF
VERBOSE

niko · 2018-05-24 08:38:11

Hi,
The problem occurs during microiterations, ok. Did you check your parameter file? Is it a standard or a modified one? In principle, the vdw radius of oxygen + the bonded terms describing H2O should be enough to avoid MM H-O dissociation.
Another possibility: how large are the ESPF charges carried by the QM atoms? Could you try your minimal example adding the Mulliken option to &espf line "External = Tinker"? It would pass Mulliken charges to Tinker instead of ESPF ones.

ArnoldBax · 2018-05-24 09:19:32

Hi,

with Mulliken charges it works fine: thanks! Indeed at the first step with ESPF charges there's a charge as large as -6, and then presumably the electrostatic attraction is enough to cause H-O dissociation. The next question would be how accurate are Mulliken charges for the QM/MM electrostatics?

I'm using standard TIP3P parameters, though I'm new to Tinker and don't know if the H-H bond is accounted for even if it is not in the connectivity of the xyz file. At any rate, I don't have this problem with water molecules only, because it also happens with a Serine MM polar hydrogen.

=============
Updated multipoles

New QM point charge value for atom 4 = -0.09700
New QM point charge value for atom 5 = 1.14200
New QM point charge value for atom 6 = -6.25900
=============

niko · 2018-05-24 09:25:47

Using any atomic charge, but the ESPF one, breaks down the unicity of the potential energy surface. In other words, atoms (either QM or MM) optimized in Molcas are moving on a different surface different from the Tinker MM one. How large is this difference is not documented and is probably system-dependent.

Regarding your ESPF charges, something is really weird with the total charge of your system which is -6! Or is it on purpose? In the molcas output, at the end of the SCF module, do you find the very same charges in the ESPF analysis part?

ArnoldBax · 2018-05-24 10:25:07

Ok, I realised Tinker was reading the wrong charges, because there was an old *.Espf.Data file in the scratch directory. Now the water dimer systems works fine, and with reasonable ESPF charges. Sorry for the mistake. Now I'm checking with the larger system.

Molcas Forum

Announcement

#1 2018-05-23 12:16:56

MM atom collapsing into QM nucleus in Tinker/Molcas

#2 2018-05-23 13:27:41

Re: MM atom collapsing into QM nucleus in Tinker/Molcas

#3 2018-05-23 14:21:49

Re: MM atom collapsing into QM nucleus in Tinker/Molcas

#4 2018-05-23 18:06:53

Re: MM atom collapsing into QM nucleus in Tinker/Molcas

#5 2018-05-24 08:38:11

Re: MM atom collapsing into QM nucleus in Tinker/Molcas

#6 2018-05-24 09:19:32

Re: MM atom collapsing into QM nucleus in Tinker/Molcas

#7 2018-05-24 09:25:47

Re: MM atom collapsing into QM nucleus in Tinker/Molcas

#8 2018-05-24 10:25:07

Re: MM atom collapsing into QM nucleus in Tinker/Molcas

Board footer