2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/gmxlib/network.h"
51 #include "gromacs/gmxlib/nrnb.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdlib/force.h"
55 #include "gromacs/mdlib/ns.h"
56 #include "gromacs/mdlib/qmmm.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/smalloc.h"
62 /* mopac interface routines */
64 F77_FUNC(domldt
, DOMLDT
) (int *nrqmat
, int labels
[], char keywords
[]);
67 F77_FUNC(domop
, DOMOP
) (int *nrqmat
, double qmcrd
[], int *nrmmat
,
68 double mmchrg
[], double mmcrd
[], double qmgrad
[],
69 double mmgrad
[], double *energy
, double qmcharges
[]);
73 void init_mopac(t_QMrec
*qm
)
75 /* initializes the mopac routines ans sets up the semiempirical
76 * computation by calling moldat(). The inline mopac routines can
77 * only perform gradient operations. If one would like to optimize a
78 * structure or find a transition state at PM3 level, gaussian is
86 if (!qm
->bSH
) /* if rerun then grad should not be done! */
88 sprintf(keywords
, "PRECISE GEO-OK CHARGE=%d GRAD MMOK ANALYT %s\n",
90 eQMmethod_names
[qm
->QMmethod
]);
94 sprintf(keywords
, "PRECISE GEO-OK CHARGE=%d SINGLET GRAD %s C.I.=(%d,%d) root=2 MECI \n",
96 eQMmethod_names
[qm
->QMmethod
],
97 qm
->CASorbitals
, qm
->CASelectrons
/2);
99 F77_FUNC(domldt
, DOMLDT
) (&qm
->nrQMatoms
, qm
->atomicnumberQM
, keywords
);
100 fprintf(stderr
, "keywords are: %s\n", keywords
);
105 real
call_mopac(t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[])
107 /* do the actual QMMM calculation using directly linked mopac subroutines
109 double /* always double as the MOPAC routines are always compiled in
110 double precission! */
111 *qmcrd
= NULL
, *qmchrg
= NULL
, *mmcrd
= NULL
, *mmchrg
= NULL
,
112 *qmgrad
, *mmgrad
= NULL
, energy
;
117 snew(qmcrd
, 3*(qm
->nrQMatoms
));
118 snew(qmgrad
, 3*(qm
->nrQMatoms
));
119 /* copy the data from qr into the arrays that are going to be used
120 * in the fortran routines of MOPAC
122 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
124 for (j
= 0; j
< DIM
; j
++)
126 qmcrd
[3*i
+j
] = (double)qm
->xQM
[i
][j
]*10;
131 /* later we will add the point charges here. There are some
132 * conceptual problems with semi-empirical QM in combination with
133 * point charges that we need to solve first....
135 gmx_fatal(FARGS
, "At present only ONIOM is allowed in combination"
136 " with MOPAC QM subroutines\n");
140 /* now compute the energy and the gradients.
143 snew(qmchrg
, qm
->nrQMatoms
);
144 F77_FUNC(domop
, DOMOP
) (&qm
->nrQMatoms
, qmcrd
, &mm
->nrMMatoms
,
145 mmchrg
, mmcrd
, qmgrad
, mmgrad
, &energy
, qmchrg
);
146 /* add the gradients to the f[] array, and also to the fshift[].
147 * the mopac gradients are in kCal/angstrom.
149 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
151 for (j
= 0; j
< DIM
; j
++)
153 f
[i
][j
] = (real
)10*CAL2JOULE
*qmgrad
[3*i
+j
];
154 fshift
[i
][j
] = (real
)10*CAL2JOULE
*qmgrad
[3*i
+j
];
157 QMener
= (real
)CAL2JOULE
*energy
;
158 /* do we do something with the mulliken charges?? */
167 real
call_mopac_SH(t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[])
169 /* do the actual SH QMMM calculation using directly linked mopac
172 double /* always double as the MOPAC routines are always compiled in
173 double precission! */
174 *qmcrd
= NULL
, *qmchrg
= NULL
, *mmcrd
= NULL
, *mmchrg
= NULL
,
175 *qmgrad
, *mmgrad
= NULL
, energy
;
181 snew(qmcrd
, 3*(qm
->nrQMatoms
));
182 snew(qmgrad
, 3*(qm
->nrQMatoms
));
183 /* copy the data from qr into the arrays that are going to be used
184 * in the fortran routines of MOPAC
186 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
188 for (j
= 0; j
< DIM
; j
++)
190 qmcrd
[3*i
+j
] = (double)qm
->xQM
[i
][j
]*10;
195 /* later we will add the point charges here. There are some
196 * conceptual problems with semi-empirical QM in combination with
197 * point charges that we need to solve first....
199 gmx_fatal(FARGS
, "At present only ONIOM is allowed in combination with MOPAC\n");
203 /* now compute the energy and the gradients.
205 snew(qmchrg
, qm
->nrQMatoms
);
207 F77_FUNC(domop
, DOMOP
) (&qm
->nrQMatoms
, qmcrd
, &mm
->nrMMatoms
,
208 mmchrg
, mmcrd
, qmgrad
, mmgrad
, &energy
, qmchrg
);
209 /* add the gradients to the f[] array, and also to the fshift[].
210 * the mopac gradients are in kCal/angstrom.
212 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
214 for (j
= 0; j
< DIM
; j
++)
216 f
[i
][j
] = (real
)10*CAL2JOULE
*qmgrad
[3*i
+j
];
217 fshift
[i
][j
] = (real
)10*CAL2JOULE
*qmgrad
[3*i
+j
];
220 QMener
= (real
)CAL2JOULE
*energy
;
225 } /* call_mopac_SH */
229 gmx_qmmm_mopac_empty
;