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,2016, 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.
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/gmxlib/network.h"
50 #include "gromacs/gmxlib/nrnb.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdlib/force.h"
54 #include "gromacs/mdlib/ns.h"
55 #include "gromacs/mdlib/qmmm.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/smalloc.h"
61 /* TODO: this should be made thread-safe */
63 /* Gaussian interface routines */
65 void init_gaussian(t_QMrec
*qm
)
69 basissets
[eQMbasisNR
] = {{0, 3, 0},
70 {0, 3, 0}, /*added for double sto-3g entry in names.c*/
84 /* using the ivec above to convert the basis read form the mdp file
85 * in a human readable format into some numbers for the gaussian
86 * route. This is necessary as we are using non standard routes to
90 /* per layer we make a new subdir for integral file, checkpoint
91 * files and such. These dirs are stored in the QMrec for
96 if (!qm
->nQMcpus
) /* this we do only once per layer
97 * as we call g01 externally
101 for (i
= 0; i
< DIM
; i
++)
103 qm
->SHbasis
[i
] = basissets
[qm
->QMbasis
][i
];
106 /* init gradually switching on of the SA */
108 /* we read the number of cpus and environment from the environment
111 buf
= getenv("GMX_QM_GAUSSIAN_NCPUS");
114 sscanf(buf
, "%d", &qm
->nQMcpus
);
120 fprintf(stderr
, "number of CPUs for gaussian = %d\n", qm
->nQMcpus
);
121 buf
= getenv("GMX_QM_GAUSSIAN_MEMORY");
124 sscanf(buf
, "%d", &qm
->QMmem
);
128 qm
->QMmem
= 50000000;
130 fprintf(stderr
, "memory for gaussian = %d\n", qm
->QMmem
);
131 buf
= getenv("GMX_QM_ACCURACY");
134 sscanf(buf
, "%d", &qm
->accuracy
);
140 fprintf(stderr
, "accuracy in l510 = %d\n", qm
->accuracy
);
142 buf
= getenv("GMX_QM_CPMCSCF");
145 sscanf(buf
, "%d", &i
);
146 qm
->cpmcscf
= (i
!= 0);
154 fprintf(stderr
, "using cp-mcscf in l1003\n");
158 fprintf(stderr
, "NOT using cp-mcscf in l1003\n");
160 buf
= getenv("GMX_QM_SA_STEP");
163 sscanf(buf
, "%d", &qm
->SAstep
);
167 /* init gradually switching on of the SA */
170 /* we read the number of cpus and environment from the environment
173 fprintf(stderr
, "Level of SA at start = %d\n", qm
->SAstep
);
174 /* punch the LJ C6 and C12 coefficients to be picked up by
175 * gaussian and usd to compute the LJ interaction between the
178 if (qm
->bTS
|| qm
->bOPT
)
180 out
= fopen("LJ.dat", "w");
181 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
183 fprintf(out
, "%3d %10.7f %10.7f\n",
184 qm
->atomicnumberQM
[i
], qm
->c6
[i
], qm
->c12
[i
]);
188 /* gaussian settings on the system */
189 buf
= getenv("GMX_QM_GAUSS_DIR");
193 qm
->gauss_dir
= gmx_strdup(buf
);
197 gmx_fatal(FARGS
, "no $GMX_QM_GAUSS_DIR, check gaussian manual\n");
200 buf
= getenv("GMX_QM_GAUSS_EXE");
203 qm
->gauss_exe
= gmx_strdup(buf
);
207 gmx_fatal(FARGS
, "no $GMX_QM_GAUSS_EXE set, check gaussian manual\n");
209 buf
= getenv("GMX_QM_MODIFIED_LINKS_DIR");
212 qm
->devel_dir
= gmx_strdup (buf
);
216 gmx_fatal(FARGS
, "no $GMX_QM_MODIFIED_LINKS_DIR, this is were the modified links reside.\n");
220 /* reactionfield, file is needed using gaussian */
221 /* rffile=fopen("rf.dat","w");*/
222 /* fprintf(rffile,"%f %f\n",fr->epsilon_r,fr->rcoulomb/BOHR2NM);*/
226 fprintf(stderr
, "gaussian initialised...\n");
231 void write_gaussian_SH_input(int step
, gmx_bool swap
,
232 t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
)
243 bSA
= (qm
->SAstep
> 0);
245 out
= fopen("input.com", "w");
246 /* write the route */
247 fprintf(out
, "%s", "%scr=input\n");
248 fprintf(out
, "%s", "%rwf=input\n");
249 fprintf(out
, "%s", "%int=input\n");
250 fprintf(out
, "%s", "%d2e=input\n");
252 * fprintf(out,"%s","%nosave\n");
254 fprintf(out
, "%s", "%chk=input\n");
255 fprintf(out
, "%s%d\n", "%mem=", qm
->QMmem
);
256 fprintf(out
, "%s%3d\n", "%nprocshare=", qm
->nQMcpus
);
258 /* use the versions of
259 * l301 that computes the interaction between MM and QM atoms.
260 * l510 that can punch the CI coefficients
261 * l701 that can do gradients on MM atoms
265 fprintf(out
, "%s%s%s",
269 fprintf(out
, "%s%s%s",
273 fprintf(out
, "%s%s%s",
278 fprintf(out
, "%s%s%s",
282 fprintf(out
, "%s%s%s",
286 /* print the nonstandard route
289 "#P nonstd\n 1/18=10,20=1,38=1/1;\n");
291 " 2/9=110,15=1,17=6,18=5,40=1/2;\n");
295 " 3/5=%d,6=%d,7=%d,25=1,32=1,43=1,94=-2/1,2,3;\n",
298 qm
->SHbasis
[2]); /*basisset stuff */
303 " 3/5=%d,6=%d,7=%d,25=1,32=1,43=0,94=-2/1,2,3;\n",
306 qm
->SHbasis
[2]); /*basisset stuff */
309 if (step
+1) /* fetch initial guess from check point file */
310 { /* hack, to alyays read from chk file!!!!! */
311 fprintf(out
, "%s%d,%s%d%s", " 4/5=1,7=6,17=",
313 "18=", qm
->CASorbitals
, "/1,5;\n");
315 else /* generate the first checkpoint file */
317 fprintf(out
, "%s%d,%s%d%s", " 4/5=0,7=6,17=",
319 "18=", qm
->CASorbitals
, "/1,5;\n");
321 /* the rest of the input depends on where the system is on the PES
323 if (swap
&& bSA
) /* make a slide to the other surface */
325 if (qm
->CASorbitals
> 6) /* use direct and no full diag */
327 fprintf(out
, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6,97=100/10;\n");
333 fprintf(out
, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6,97=100/10;\n",
335 if (mm
->nrMMatoms
> 0)
337 fprintf(out
, " 7/7=1,16=-2,30=1/1;\n");
339 fprintf(out
, " 11/31=1,42=1,45=1/1;\n");
340 fprintf(out
, " 10/6=1,10=700006,28=2,29=1,31=1,97=100/3;\n");
341 fprintf(out
, " 7/30=1/16;\n 99/10=4/99;\n");
345 fprintf(out
, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6,97=100/10;\n",
347 fprintf(out
, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
351 else if (bSA
) /* do a "state-averaged" CAS calculation */
353 if (qm
->CASorbitals
> 6) /* no full diag */
355 fprintf(out
, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6/10;\n");
361 fprintf(out
, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6/10;\n",
363 if (mm
->nrMMatoms
> 0)
365 fprintf(out
, " 7/7=1,16=-2,30=1/1;\n");
367 fprintf(out
, " 11/31=1,42=1,45=1/1;\n");
368 fprintf(out
, " 10/6=1,10=700006,28=2,29=1,31=1/3;\n");
369 fprintf(out
, " 7/30=1/16;\n 99/10=4/99;\n");
373 fprintf(out
, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6/10;\n",
375 fprintf(out
, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
379 else if (swap
) /* do a "swapped" CAS calculation */
381 if (qm
->CASorbitals
> 6)
383 fprintf(out
, " 5/5=2,16=-2,17=0,28=2,32=2,38=6,97=100/10;\n");
387 fprintf(out
, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6,97=100/10;\n",
390 fprintf(out
, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
392 else /* do a "normal" CAS calculation */
394 if (qm
->CASorbitals
> 6)
396 fprintf(out
, " 5/5=2,16=-2,17=0,28=2,32=2,38=6/10;\n");
400 fprintf(out
, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6/10;\n",
403 fprintf(out
, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
405 fprintf(out
, "\ninput-file generated by gromacs\n\n");
406 fprintf(out
, "%2d%2d\n", qm
->QMcharge
, qm
->multiplicity
);
407 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
409 fprintf(out
, "%3d %10.7f %10.7f %10.7f\n",
410 qm
->atomicnumberQM
[i
],
411 qm
->xQM
[i
][XX
]/BOHR2NM
,
412 qm
->xQM
[i
][YY
]/BOHR2NM
,
413 qm
->xQM
[i
][ZZ
]/BOHR2NM
);
415 /* MM point charge data */
416 if (QMMMrec
->QMMMscheme
!= eQMMMschemeoniom
&& mm
->nrMMatoms
)
419 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
421 fprintf(out
, "%10.7f %10.7f %10.7f %8.4f\n",
422 mm
->xMM
[i
][XX
]/BOHR2NM
,
423 mm
->xMM
[i
][YY
]/BOHR2NM
,
424 mm
->xMM
[i
][ZZ
]/BOHR2NM
,
428 if (bSA
) /* put the SA coefficients at the end of the file */
430 fprintf(out
, "\n%10.8f %10.8f\n",
431 qm
->SAstep
*0.5/qm
->SAsteps
,
432 1-qm
->SAstep
*0.5/qm
->SAsteps
);
433 fprintf(stderr
, "State Averaging level = %d/%d\n", qm
->SAstep
, qm
->SAsteps
);
437 } /* write_gaussian_SH_input */
439 void write_gaussian_input(int step
, t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
)
449 out
= fopen("input.com", "w");
450 /* write the route */
452 if (qm
->QMmethod
>= eQMmethodRHF
)
464 fprintf(out
, "%s%3d\n",
465 "%nprocshare=", qm
->nQMcpus
);
467 fprintf(out
, "%s%d\n",
469 /* use the modified links that include the LJ contribution at the QM level */
470 if (qm
->bTS
|| qm
->bOPT
)
472 fprintf(out
, "%s%s%s",
473 "%subst l701 ", qm
->devel_dir
, "/l701_LJ\n");
474 fprintf(out
, "%s%s%s",
475 "%subst l301 ", qm
->devel_dir
, "/l301_LJ\n");
479 fprintf(out
, "%s%s%s",
480 "%subst l701 ", qm
->devel_dir
, "/l701\n");
481 fprintf(out
, "%s%s%s",
482 "%subst l301 ", qm
->devel_dir
, "/l301\n");
484 fprintf(out
, "%s%s%s",
485 "%subst l9999 ", qm
->devel_dir
, "/l9999\n");
496 if (qm
->QMmethod
== eQMmethodB3LYPLAN
)
499 "B3LYP/GEN Pseudo=Read");
504 eQMmethod_names
[qm
->QMmethod
]);
506 if (qm
->QMmethod
>= eQMmethodRHF
)
508 if (qm
->QMmethod
== eQMmethodCASSCF
)
510 /* in case of cas, how many electrons and orbitals do we need?
512 fprintf(out
, "(%d,%d)",
513 qm
->CASelectrons
, qm
->CASorbitals
);
516 eQMbasis_names
[qm
->QMbasis
]);
519 if (QMMMrec
->QMMMscheme
== eQMMMschemenormal
&& mm
->nrMMatoms
)
524 if (step
|| qm
->QMmethod
== eQMmethodCASSCF
)
526 /* fetch guess from checkpoint file, always for CASSCF */
527 fprintf(out
, "%s", " guess=read");
529 fprintf(out
, "\nNosymm units=bohr\n");
533 fprintf(out
, "OPT=(Redundant,TS,noeigentest,ModRedundant) Punch=(Coord,Derivatives) ");
537 fprintf(out
, "OPT=(Redundant,ModRedundant) Punch=(Coord,Derivatives) ");
541 fprintf(out
, "FORCE Punch=(Derivatives) ");
543 fprintf(out
, "iop(3/33=1)\n\n");
544 fprintf(out
, "input-file generated by gromacs\n\n");
545 fprintf(out
, "%2d%2d\n", qm
->QMcharge
, qm
->multiplicity
);
546 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
548 fprintf(out
, "%3d %10.7f %10.7f %10.7f\n",
549 qm
->atomicnumberQM
[i
],
550 qm
->xQM
[i
][XX
]/BOHR2NM
,
551 qm
->xQM
[i
][YY
]/BOHR2NM
,
552 qm
->xQM
[i
][ZZ
]/BOHR2NM
);
555 /* Pseudo Potential and ECP are included here if selected (MEthod suffix LAN) */
556 if (qm
->QMmethod
== eQMmethodB3LYPLAN
)
559 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
561 if (qm
->atomicnumberQM
[i
] < 21)
563 fprintf(out
, "%d ", i
+1);
566 fprintf(out
, "\n%s\n****\n", eQMbasis_names
[qm
->QMbasis
]);
568 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
570 if (qm
->atomicnumberQM
[i
] > 21)
572 fprintf(out
, "%d ", i
+1);
575 fprintf(out
, "\n%s\n****\n\n", "lanl2dz");
577 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
579 if (qm
->atomicnumberQM
[i
] > 21)
581 fprintf(out
, "%d ", i
+1);
584 fprintf(out
, "\n%s\n", "lanl2dz");
589 /* MM point charge data */
590 if (QMMMrec
->QMMMscheme
!= eQMMMschemeoniom
&& mm
->nrMMatoms
)
592 fprintf(stderr
, "nr mm atoms in gaussian.c = %d\n", mm
->nrMMatoms
);
594 if (qm
->bTS
|| qm
->bOPT
)
596 /* freeze the frontier QM atoms and Link atoms. This is
597 * important only if a full QM subsystem optimization is done
598 * with a frozen MM environmeent. For dynamics, or gromacs's own
599 * optimization routines this is not important.
601 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
603 if (qm
->frontatoms
[i
])
605 fprintf(out
, "%d F\n", i
+1); /* counting from 1 */
608 /* MM point charges include LJ parameters in case of QM optimization
610 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
612 fprintf(out
, "%10.7f %10.7f %10.7f %8.4f 0.0 %10.7f %10.7f\n",
613 mm
->xMM
[i
][XX
]/BOHR2NM
,
614 mm
->xMM
[i
][YY
]/BOHR2NM
,
615 mm
->xMM
[i
][ZZ
]/BOHR2NM
,
617 mm
->c6
[i
], mm
->c12
[i
]);
623 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
625 fprintf(out
, "%10.7f %10.7f %10.7f %8.4f\n",
626 mm
->xMM
[i
][XX
]/BOHR2NM
,
627 mm
->xMM
[i
][YY
]/BOHR2NM
,
628 mm
->xMM
[i
][ZZ
]/BOHR2NM
,
638 } /* write_gaussian_input */
640 real
read_gaussian_output(rvec QMgrad
[], rvec MMgrad
[], t_QMrec
*qm
, t_MMrec
*mm
)
651 in
= fopen("fort.7", "r");
655 /* in case of an optimization, the coordinates are printed in the
656 * fort.7 file first, followed by the energy, coordinates and (if
657 * required) the CI eigenvectors.
659 if (qm
->bTS
|| qm
->bOPT
)
661 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
663 if (NULL
== fgets(buf
, 300, in
))
665 gmx_fatal(FARGS
, "Error reading Gaussian output - not enough atom lines?");
669 sscanf(buf
, "%d %lf %lf %lf\n",
675 sscanf(buf
, "%d %f %f %f\n",
681 for (j
= 0; j
< DIM
; j
++)
683 qm
->xQM
[i
][j
] *= BOHR2NM
;
687 /* the next line is the energy and in the case of CAS, the energy
688 * difference between the two states.
690 if (NULL
== fgets(buf
, 300, in
))
692 gmx_fatal(FARGS
, "Error reading Gaussian output");
696 sscanf(buf
, "%lf\n", &QMener
);
698 sscanf(buf
, "%f\n", &QMener
);
700 /* next lines contain the gradients of the QM atoms */
701 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
703 if (NULL
== fgets(buf
, 300, in
))
705 gmx_fatal(FARGS
, "Error reading Gaussian output");
708 sscanf(buf
, "%lf %lf %lf\n",
713 sscanf(buf
, "%f %f %f\n",
719 /* the next lines are the gradients of the MM atoms */
720 if (qm
->QMmethod
>= eQMmethodRHF
)
722 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
724 if (NULL
== fgets(buf
, 300, in
))
726 gmx_fatal(FARGS
, "Error reading Gaussian output");
729 sscanf(buf
, "%lf %lf %lf\n",
734 sscanf(buf
, "%f %f %f\n",
745 real
read_gaussian_SH_output(rvec QMgrad
[], rvec MMgrad
[], int step
, t_QMrec
*qm
, t_MMrec
*mm
)
756 in
= fopen("fort.7", "r");
757 /* first line is the energy and in the case of CAS, the energy
758 * difference between the two states.
760 if (NULL
== fgets(buf
, 300, in
))
762 gmx_fatal(FARGS
, "Error reading Gaussian output");
766 sscanf(buf
, "%lf %lf\n", &QMener
, &DeltaE
);
768 sscanf(buf
, "%f %f\n", &QMener
, &DeltaE
);
771 /* switch on/off the State Averaging */
773 if (DeltaE
> qm
->SAoff
)
780 else if (DeltaE
< qm
->SAon
|| (qm
->SAstep
> 0))
782 if (qm
->SAstep
< qm
->SAsteps
)
789 fprintf(stderr
, "Gap = %5f,SA = %3d\n", DeltaE
, (qm
->SAstep
> 0));
790 /* next lines contain the gradients of the QM atoms */
791 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
793 if (NULL
== fgets(buf
, 300, in
))
795 gmx_fatal(FARGS
, "Error reading Gaussian output");
799 sscanf(buf
, "%lf %lf %lf\n",
804 sscanf(buf
, "%f %f %f\n",
810 /* the next lines, are the gradients of the MM atoms */
812 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
814 if (NULL
== fgets(buf
, 300, in
))
816 gmx_fatal(FARGS
, "Error reading Gaussian output");
819 sscanf(buf
, "%lf %lf %lf\n",
824 sscanf(buf
, "%f %f %f\n",
831 /* the next line contains the two CI eigenvector elements */
832 if (NULL
== fgets(buf
, 300, in
))
834 gmx_fatal(FARGS
, "Error reading Gaussian output");
838 sscanf(buf
, "%d", &qm
->CIdim
);
839 snew(qm
->CIvec1
, qm
->CIdim
);
840 snew(qm
->CIvec1old
, qm
->CIdim
);
841 snew(qm
->CIvec2
, qm
->CIdim
);
842 snew(qm
->CIvec2old
, qm
->CIdim
);
846 /* before reading in the new current CI vectors, copy the current
847 * CI vector into the old one.
849 for (i
= 0; i
< qm
->CIdim
; i
++)
851 qm
->CIvec1old
[i
] = qm
->CIvec1
[i
];
852 qm
->CIvec2old
[i
] = qm
->CIvec2
[i
];
856 for (i
= 0; i
< qm
->CIdim
; i
++)
858 if (NULL
== fgets(buf
, 300, in
))
860 gmx_fatal(FARGS
, "Error reading Gaussian output");
863 sscanf(buf
, "%lf\n", &qm
->CIvec1
[i
]);
865 sscanf(buf
, "%f\n", &qm
->CIvec1
[i
]);
869 for (i
= 0; i
< qm
->CIdim
; i
++)
871 if (NULL
== fgets(buf
, 300, in
))
873 gmx_fatal(FARGS
, "Error reading Gaussian output");
876 sscanf(buf
, "%lf\n", &qm
->CIvec2
[i
]);
878 sscanf(buf
, "%f\n", &qm
->CIvec2
[i
]);
885 real
inproduct(real
*a
, real
*b
, int n
)
892 /* computes the inner product between two vectors (a.b), both of
893 * which have length n.
895 for (i
= 0; i
< n
; i
++)
902 int hop(int step
, t_QMrec
*qm
)
907 d11
= 0.0, d12
= 0.0, d21
= 0.0, d22
= 0.0;
909 /* calculates the inproduct between the current Ci vector and the
910 * previous CI vector. A diabatic hop will be made if d12 and d21
911 * are much bigger than d11 and d22. In that case hop returns true,
912 * otherwise it returns false.
914 if (step
) /* only go on if more than one step has been done */
916 d11
= inproduct(qm
->CIvec1
, qm
->CIvec1old
, qm
->CIdim
);
917 d12
= inproduct(qm
->CIvec1
, qm
->CIvec2old
, qm
->CIdim
);
918 d21
= inproduct(qm
->CIvec2
, qm
->CIvec1old
, qm
->CIdim
);
919 d22
= inproduct(qm
->CIvec2
, qm
->CIvec2old
, qm
->CIdim
);
921 fprintf(stderr
, "-------------------\n");
922 fprintf(stderr
, "d11 = %13.8f\n", d11
);
923 fprintf(stderr
, "d12 = %13.8f\n", d12
);
924 fprintf(stderr
, "d21 = %13.8f\n", d21
);
925 fprintf(stderr
, "d22 = %13.8f\n", d22
);
926 fprintf(stderr
, "-------------------\n");
928 if ((fabs(d12
) > 0.5) && (fabs(d21
) > 0.5))
936 void do_gaussian(int step
, char *exe
)
941 /* make the call to the gaussian binary through system()
942 * The location of the binary will be picked up from the
943 * environment using getenv().
945 if (step
) /* hack to prevent long inputfiles */
947 sprintf(buf
, "%s < %s > %s",
954 sprintf(buf
, "%s < %s > %s",
959 fprintf(stderr
, "Calling '%s'\n", buf
);
960 if (system(buf
) != 0)
962 gmx_fatal(FARGS
, "Call to '%s' failed\n", buf
);
966 real
call_gaussian(t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[])
968 /* normal gaussian jobs */
981 sprintf(exe
, "%s/%s", qm
->gauss_dir
, qm
->gauss_exe
);
982 snew(QMgrad
, qm
->nrQMatoms
);
983 snew(MMgrad
, mm
->nrMMatoms
);
985 write_gaussian_input(step
, fr
, qm
, mm
);
986 do_gaussian(step
, exe
);
987 QMener
= read_gaussian_output(QMgrad
, MMgrad
, qm
, mm
);
988 /* put the QMMM forces in the force array and to the fshift
990 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
992 for (j
= 0; j
< DIM
; j
++)
994 f
[i
][j
] = HARTREE_BOHR2MD
*QMgrad
[i
][j
];
995 fshift
[i
][j
] = HARTREE_BOHR2MD
*QMgrad
[i
][j
];
998 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
1000 for (j
= 0; j
< DIM
; j
++)
1002 f
[i
+qm
->nrQMatoms
][j
] = HARTREE_BOHR2MD
*MMgrad
[i
][j
];
1003 fshift
[i
+qm
->nrQMatoms
][j
] = HARTREE_BOHR2MD
*MMgrad
[i
][j
];
1006 QMener
= QMener
*HARTREE2KJ
*AVOGADRO
;
1011 } /* call_gaussian */
1013 real
call_gaussian_SH(t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[])
1015 /* a gaussian call routine intended for doing diabatic surface
1016 * "sliding". See the manual for the theoretical background of this
1026 swapped
= FALSE
; /* handle for identifying the current PES */
1028 swap
= FALSE
; /* the actual swap */
1037 sprintf(exe
, "%s/%s", qm
->gauss_dir
, qm
->gauss_exe
);
1038 /* hack to do ground state simulations */
1042 buf
= getenv("GMX_QM_GROUND_STATE");
1045 sscanf(buf
, "%d", &state
);
1059 /* copy the QMMMrec pointer */
1060 snew(QMgrad
, qm
->nrQMatoms
);
1061 snew(MMgrad
, mm
->nrMMatoms
);
1062 /* at step 0 there should be no SA */
1065 /* temporray set to step + 1, since there is a chk start */
1066 write_gaussian_SH_input(step
, swapped
, fr
, qm
, mm
);
1068 do_gaussian(step
, exe
);
1069 QMener
= read_gaussian_SH_output(QMgrad
, MMgrad
, step
, qm
, mm
);
1071 /* check for a surface hop. Only possible if we were already state
1078 swap
= (step
&& hop(step
, qm
));
1081 else /* already on the other surface, so check if we go back */
1083 swap
= (step
&& hop(step
, qm
));
1084 swapped
= !swap
; /* so swapped shoud be false again */
1086 if (swap
) /* change surface, so do another call */
1088 write_gaussian_SH_input(step
, swapped
, fr
, qm
, mm
);
1089 do_gaussian(step
, exe
);
1090 QMener
= read_gaussian_SH_output(QMgrad
, MMgrad
, step
, qm
, mm
);
1093 /* add the QMMM forces to the gmx force array and fshift
1095 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1097 for (j
= 0; j
< DIM
; j
++)
1099 f
[i
][j
] = HARTREE_BOHR2MD
*QMgrad
[i
][j
];
1100 fshift
[i
][j
] = HARTREE_BOHR2MD
*QMgrad
[i
][j
];
1103 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
1105 for (j
= 0; j
< DIM
; j
++)
1107 f
[i
+qm
->nrQMatoms
][j
] = HARTREE_BOHR2MD
*MMgrad
[i
][j
];
1108 fshift
[i
+qm
->nrQMatoms
][j
] = HARTREE_BOHR2MD
*MMgrad
[i
][j
];
1111 QMener
= QMener
*HARTREE2KJ
*AVOGADRO
;
1112 fprintf(stderr
, "step %5d, SA = %5d, swap = %5d\n",
1113 step
, (qm
->SAstep
> 0), swapped
);
1118 } /* call_gaussian_SH */
1120 /* end of gaussian sub routines */
1124 gmx_qmmm_gaussian_empty
;