[gmx-users] QM/MM simulation: point-charge, pbc and bOPT=yes problems

[Bence Hégely]

Dear Gromacs users!
I'm trying to interface gromacs (5.0) with orca (2.9.1), although i
encountered some problems. I'm testing a tripeptide system which i found
here: http://www.picb.ac.cn/~liwenjin/QMMM_simulations/index.html , but i'm
using a different .mdp file, with these parameters (i copied the only
working configuration):
integrator          = steep
dt                  =  0.002    ; ps !
nsteps              =   10000
nstlist             =  10
ns_type             =  grid
rlist               =  3.0
rcoulomb            =  3.0
coulombtype         = cut-off
vdwtype             = cut-off
rvdw                = 3.0
pbc                 = no
periodic_molecules  =  no
constraints         = none
energygrps          = QMatoms MMatoms
cutoff-scheme       = group

; QM/MM calculation stuff
QMMM = yes
QMMM-grps = QMatoms
QMmethod = rhf
QMbasis = 3-21G
QMMMscheme = normal
QMcharge = 0
QMmult = 1
bOPT                     = no
bTS                      = no
SH                       = no

;
;       Energy minimizing stuff
;
emtol               =  60   ; minimization thresold (kj/mol.nm-1)    1
hartree/bohr= 49614.75241 kj/mol.nm-1  1 kj/mol.nm-1=2.01553e-5 hartree/bohr
emstep              =  0.01  ; minimization step in nm

with the .ORCAINFO file of
! LDA cc-(p)VDZ

%LJcoefficients "peptide.LJ"

%pointcharges "peptide.pc"

%geom

end

After i modified the ffnonbonded.itp and atomtypes.atp for the dummy atom,
rebuilt the gromacs and the grompp and mdrun commands were executed with
the following:
grompp -p peptide.top -c peptide.gro -f peptide.mdp -n peptide.ndx -o
peptide.tpr
mdrun -s peptide.tpr -c peptide.gro -o peptide.trr -e peptide.edr -g
peptide.log

The first thing is, that the default neighbor searching algorithm -
cutoff-scheme=verlet - isn't producing any point charges for ORCA and the
job ends with a segmentation fault, although the group cutoff scheme works
well. It is not clear to me why this is happening, but i suspect the
problem lies in that the qm/mm subrutines are quite old (~2003?), and the
verlet scheme introduced much later in gromacs 4.6 (2013).

The second problem that i have encountered is with the pbc options: if i'm
using pbc=xyz the job terminates with a segmentation fault (core dumped)
after orca terminates normally. After i looked for some explanation about
coordinate updating on the qmmm.c source file I saw some comments that the
update_QMMMrec function not working properly without pbc, even though it
works just fine with pbc=no. Little bit confused of that, although the
gromacs was tested on the regression package and all the tests passed.

The third problem is with the bOPT = yes option. When i pass the
optimization to the orca, after the first optimization cycle (before orca
terminates normally)  the tripeptide's qm region forms an unrealistic
geometry. At the input of the second orca run, the qm region scatters, the
atoms move apart around 10 angströms from each other and gives the
following error message:
Error (ORCA_GSTEP): The lambda equations have not converged

although, the job continues and at the third orca run gives the following:
Calling '/export/home/hegely/Programz/orca_2_9_1_linux_x86-64/orca
peptide.inp >> peptide.out'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!                        FATAL ERROR ENCOUNTERED               !!!
!!!                        -----------------------               !!!
!!!                          I/O OPERATION FAILED                !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ABORTING THE RUN

and the job runs onward, giving the same error message over and over until
i kill it. In the gromacs .log file no error messages are presented and i
don't really have a clue of what's going on. Maybe the point charges
destroy the geometry of the peptide?

I think i checked all the qm/mm, orca searching results in the gmx-users
list, but didn't find the solutions. Sorry if i missed something out, but i
would appreciate any help!

Bence Hégely