<div dir="ltr">Hi all, I have a iron(III) atom with its 6 coordinated water molecules (QM region) surrounded by 2141 MM water molecules. I started with a geometry optimization at the QM/MM level and that goes well (although only 1.5 charge is localised on the iron). I then start a MD run and things get messy. The error is about one QM atom (the iron) being too close to a MM atom. Then the MM energies go to the roof, I have NaN gradient and evrything crashes. I don't understand how this works. First there are the QM waters and other MM water molecules between the two "too close" atoms so I don't think it's not a "real" geometry problem. I have tried to change the size of the QM box, make the QM electrostatics periodic (initially set to none periodic), try MT and multipole for the non periodic Psolver, change the number of gaussian in GEEP_LIB, use 1 grid .. Nothing seems to work. I have notice that the partial charge on the iron atom changes quite a lot between 2 steps and when it's close to 0, the MM potential becomes huge. I don't know why this is so however, and even if my iron goes neutral, it should not fly into a MM water ... 1st step MULLIKEN POPULATION ANALYSIS # Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment 1 Fe 1 9.185313 5.416885 1.397803 3.768428 2 O 2 3.479264 3.199820 -0.679084 0.279444 3 O 2 3.454824 3.180460 -0.635284 0.274363 4 O 2 3.431022 3.274684 -0.705706 0.156338 5 O 2 3.463094 3.224418 -0.687511 0.238676 6 O 2 3.426629 3.272281 -0.698910 0.154349 7 O 2 3.412406 3.254903 -0.667308 0.157503 8 H 3 0.242603 0.240340 0.517057 0.002264 9 H 3 0.246114 0.251097 0.502789 -0.004983 10 H 3 0.251669 0.257771 0.490560 -0.006102 11 H 3 0.245008 0.247299 0.507693 -0.002291 12 H 3 0.266848 0.270467 0.462684 -0.003619 13 H 3 0.266236 0.272343 0.461420 -0.006107 14 H 3 0.258389 0.257856 0.483755 0.000534 15 H 3 0.268412 0.266475 0.465113 0.001937 16 H 3 0.271005 0.274766 0.454229 -0.003762 17 H 3 0.273854 0.272849 0.453298 0.001005 18 H 3 0.273299 0.277448 0.449253 -0.004148 19 H 3 0.284011 0.287840 0.428148 -0.003829 # Total charge and spin 33.000000 28.000000 3.000000 5.000000 [...] 3rd step MULLIKEN POPULATION ANALYSIS # Atom Element Kind Atomic population (alpha,beta) Net charge Spin moment 1 Fe 1 9.000001 6.995986 0.004012 2.004015 2 O 2 3.564147 2.829655 -0.393802 0.734492 3 O 2 3.530082 2.891696 -0.421778 0.638385 4 O 2 3.477096 3.065067 -0.542163 0.412029 5 O 2 3.511836 2.968556 -0.480392 0.543280 6 O 2 3.497889 2.993893 -0.491782 0.503997 7 O 2 3.448838 3.094323 -0.543160 0.354515 8 H 3 0.219325 0.241090 0.539585 -0.021765 9 H 3 0.218180 0.242569 0.539251 -0.024390 10 H 3 0.232233 0.250945 0.516823 -0.018712 11 H 3 0.232325 0.250821 0.516853 -0.018496 12 H 3 0.250804 0.261371 0.487824 -0.010567 13 H 3 0.238978 0.258239 0.502783 -0.019260 14 H 3 0.255165 0.272122 0.472713 -0.016956 15 H 3 0.270316 0.281142 0.448542 -0.010826 16 H 3 0.247369 0.261604 0.491027 -0.014235 17 H 3 0.254925 0.269922 0.475153 -0.014996 18 H 3 0.267309 0.277233 0.455457 -0.009924 19 H 3 0.283180 0.293767 0.423053 -0.010586 # Total charge and spin 33.000000 28.000000 3.000000 5.000000 [...] QMMM| Information on the QM/MM Electrostatic Potential: QMMM| QM/MM Coupling computed collocating the Gaussian Potential Functions. WARNING| Particles: 512 1 at distance [au]: 0.58442032 less than: 0.95400000; increase EMAX_SPLINE. ******************************************************************************** * Value of r in Input = 0.584420 not in the spline range. Using = 0.901013 * ******************************************************************************** MM DIPOLE [BERRY PHASE](A.U.)| -9.832403 4.451024 19.766879 MM DIPOLE [BERRY PHASE](Debye)| -24.991473 11.313373 50.242390 MM DIPOLE [BERRY PHASE] DERIVATIVE(A.U.)| 0.000464 -0.003523 -0.000552 ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.): 1560181.245483303442597 The most curious thing is that if I switch back from MD run to Geo_OPT run, all is good, the charge on the iron stay stable and I can run it for as long as I want. Then it would seem that my forcefield parameters are good, otherwise the MM energy would not be stable there either. I'm assuming this is is a silly mistake in my input but I can't work it out and I was hoping one of you would be able to help me out. Thank you, G </div>