3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
39 #ifdef GMX_QMMM_GAUSSIAN
64 #include "gmx_fatal.h"
69 /* TODO: this should be made thread-safe */
71 /* Gaussian interface routines */
73 void init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
76 *rffile=NULL,*out=NULL;
78 basissets[eQMbasisNR]={{0,3,0},
79 {0,3,0},/*added for double sto-3g entry in names.c*/
93 /* using the ivec above to convert the basis read form the mdp file
94 * in a human readable format into some numbers for the gaussian
95 * route. This is necessary as we are using non standard routes to
99 /* per layer we make a new subdir for integral file, checkpoint
100 * files and such. These dirs are stored in the QMrec for
105 if(!qm->nQMcpus){ /* this we do only once per layer
106 * as we call g01 externally
110 qm->SHbasis[i]=basissets[qm->QMbasis][i];
112 /* init gradually switching on of the SA */
114 /* we read the number of cpus and environment from the environment
117 buf = getenv("NCPUS");
119 sscanf(buf,"%d",&qm->nQMcpus);
122 fprintf(stderr,"number of CPUs for gaussian = %d\n",qm->nQMcpus);
125 sscanf(buf,"%d",&qm->QMmem);
128 fprintf(stderr,"memory for gaussian = %d\n",qm->QMmem);
131 sscanf(buf,"%d",&qm->accuracy);
134 fprintf(stderr,"accuracy in l510 = %d\n",qm->accuracy);
136 buf = getenv("CPMCSCF");
140 qm->cpmcscf = (i!=0);
145 fprintf(stderr,"using cp-mcscf in l1003\n");
147 fprintf(stderr,"NOT using cp-mcscf in l1003\n");
148 buf = getenv("SASTEP");
150 sscanf(buf,"%d",&qm->SAstep);
152 /* init gradually switching on of the SA */
154 /* we read the number of cpus and environment from the environment
157 fprintf(stderr,"Level of SA at start = %d\n",qm->SAstep);
158 /* punch the LJ C6 and C12 coefficients to be picked up by
159 * gaussian and usd to compute the LJ interaction between the
162 if(qm->bTS||qm->bOPT){
163 out = fopen("LJ.dat","w");
164 for(i=0;i<qm->nrQMatoms;i++){
167 fprintf(out,"%3d %10.7lf %10.7lf\n",
168 qm->atomicnumberQM[i],qm->c6[i],qm->c12[i]);
170 fprintf(out,"%3d %10.7f %10.7f\n",
171 qm->atomicnumberQM[i],qm->c6[i],qm->c12[i]);
176 /* gaussian settings on the system */
177 buf = getenv("GAUSS_DIR");
178 fprintf(stderr,"%s",buf);
181 qm->gauss_dir=strdup(buf);
184 gmx_fatal(FARGS,"no $GAUSS_DIR, check gaussian manual\n");
186 buf = getenv("GAUSS_EXE");
188 qm->gauss_exe=strdup(buf);
191 gmx_fatal(FARGS,"no $GAUSS_EXE, check gaussian manual\n");
192 buf = getenv("DEVEL_DIR");
194 qm->devel_dir = strdup (buf);
197 gmx_fatal(FARGS,"no $DEVEL_DIR, this is were the modified links reside.\n");
200 /* reactionfield, file is needed using gaussian */
201 /* rffile=fopen("rf.dat","w");*/
202 /* fprintf(rffile,"%f %f\n",fr->epsilon_r,fr->rcoulomb/BOHR2NM);*/
206 fprintf(stderr,"gaussian initialised...\n");
211 void write_gaussian_SH_input(int step,gmx_bool swap,
212 t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
223 bSA = (qm->SAstep>0);
225 out = fopen("input.com","w");
226 /* write the route */
227 fprintf(out,"%s","%scr=input\n");
228 fprintf(out,"%s","%rwf=input\n");
229 fprintf(out,"%s","%int=input\n");
230 fprintf(out,"%s","%d2e=input\n");
232 * fprintf(out,"%s","%nosave\n");
234 fprintf(out,"%s","%chk=input\n");
235 fprintf(out,"%s%d\n","%mem=",qm->QMmem);
236 fprintf(out,"%s%3d\n","%nprocshare=",qm->nQMcpus);
238 /* use the versions of
239 * l301 that computes the interaction between MM and QM atoms.
240 * l510 that can punch the CI coefficients
241 * l701 that can do gradients on MM atoms
245 fprintf(out,"%s%s%s",
249 fprintf(out,"%s%s%s",
253 fprintf(out,"%s%s%s",
258 fprintf(out,"%s%s%s",
262 fprintf(out,"%s%s%s",
266 /* print the nonstandard route
269 "#P nonstd\n 1/18=10,20=1,38=1/1;\n");
271 " 2/9=110,15=1,17=6,18=5,40=1/2;\n");
274 " 3/5=%d,6=%d,7=%d,25=1,32=1,43=1,94=-2/1,2,3;\n",
277 qm->SHbasis[2]); /*basisset stuff */
280 " 3/5=%d,6=%d,7=%d,25=1,32=1,43=0,94=-2/1,2,3;\n",
283 qm->SHbasis[2]); /*basisset stuff */
285 if (step+1) /* fetch initial guess from check point file */
286 /* hack, to alyays read from chk file!!!!! */
287 fprintf(out,"%s%d,%s%d%s"," 4/5=1,7=6,17=",
289 "18=",qm->CASorbitals,"/1,5;\n");
290 else /* generate the first checkpoint file */
291 fprintf(out,"%s%d,%s%d%s"," 4/5=0,7=6,17=",
293 "18=",qm->CASorbitals,"/1,5;\n");
294 /* the rest of the input depends on where the system is on the PES
296 if(swap && bSA){ /* make a slide to the other surface */
297 if(qm->CASorbitals>6){ /* use direct and no full diag */
298 fprintf(out," 5/5=2,16=-2,17=10000000,28=2,32=2,38=6,97=100/10;\n");
302 fprintf(out," 5/5=2,6=%d,17=31000200,28=2,32=2,38=6,97=100/10;\n",
305 fprintf(out," 7/7=1,16=-2,30=1/1;\n");
306 fprintf(out," 11/31=1,42=1,45=1/1;\n");
307 fprintf(out," 10/6=1,10=700006,28=2,29=1,31=1,97=100/3;\n");
308 fprintf(out," 7/30=1/16;\n 99/10=4/99;\n");
311 fprintf(out," 5/5=2,6=%d,17=11000000,28=2,32=2,38=6,97=100/10;\n",
313 fprintf(out," 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
317 else if(bSA){ /* do a "state-averaged" CAS calculation */
318 if(qm->CASorbitals>6){ /* no full diag */
319 fprintf(out," 5/5=2,16=-2,17=10000000,28=2,32=2,38=6/10;\n");
323 fprintf(out," 5/5=2,6=%d,17=31000200,28=2,32=2,38=6/10;\n",
326 fprintf(out," 7/7=1,16=-2,30=1/1;\n");
327 fprintf(out," 11/31=1,42=1,45=1/1;\n");
328 fprintf(out," 10/6=1,10=700006,28=2,29=1,31=1/3;\n");
329 fprintf(out," 7/30=1/16;\n 99/10=4/99;\n");
332 fprintf(out," 5/5=2,6=%d,17=11000000,28=2,32=2,38=6/10;\n",
334 fprintf(out," 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
338 else if(swap){/* do a "swapped" CAS calculation */
339 if(qm->CASorbitals>6)
340 fprintf(out," 5/5=2,16=-2,17=0,28=2,32=2,38=6,97=100/10;\n");
342 fprintf(out," 5/5=2,6=%d,17=1000000,28=2,32=2,38=6,97=100/10;\n",
344 fprintf(out," 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
346 else {/* do a "normal" CAS calculation */
347 if(qm->CASorbitals>6)
348 fprintf(out," 5/5=2,16=-2,17=0,28=2,32=2,38=6/10;\n");
350 fprintf(out," 5/5=2,6=%d,17=1000000,28=2,32=2,38=6/10;\n",
352 fprintf(out," 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
354 fprintf(out, "\ninput-file generated by gromacs\n\n");
355 fprintf(out,"%2d%2d\n",qm->QMcharge,qm->multiplicity);
356 for (i=0;i<qm->nrQMatoms;i++){
358 fprintf(out,"%3d %10.7lf %10.7lf %10.7lf\n",
359 qm->atomicnumberQM[i],
360 qm->xQM[i][XX]/BOHR2NM,
361 qm->xQM[i][YY]/BOHR2NM,
362 qm->xQM[i][ZZ]/BOHR2NM);
364 fprintf(out,"%3d %10.7f %10.7f %10.7f\n",
365 qm->atomicnumberQM[i],
366 qm->xQM[i][XX]/BOHR2NM,
367 qm->xQM[i][YY]/BOHR2NM,
368 qm->xQM[i][ZZ]/BOHR2NM);
371 /* MM point charge data */
372 if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
374 for(i=0;i<mm->nrMMatoms;i++){
376 fprintf(out,"%10.7lf %10.7lf %10.7lf %8.4lf\n",
377 mm->xMM[i][XX]/BOHR2NM,
378 mm->xMM[i][YY]/BOHR2NM,
379 mm->xMM[i][ZZ]/BOHR2NM,
382 fprintf(out,"%10.7f %10.7f %10.7f %8.4f\n",
383 mm->xMM[i][XX]/BOHR2NM,
384 mm->xMM[i][YY]/BOHR2NM,
385 mm->xMM[i][ZZ]/BOHR2NM,
390 if(bSA) {/* put the SA coefficients at the end of the file */
392 fprintf(out,"\n%10.8lf %10.8lf\n",
393 qm->SAstep*0.5/qm->SAsteps,
394 1-qm->SAstep*0.5/qm->SAsteps);
396 fprintf(out,"\n%10.8f %10.8f\n",
397 qm->SAstep*0.5/qm->SAsteps,
398 1-qm->SAstep*0.5/qm->SAsteps);
400 fprintf(stderr,"State Averaging level = %d/%d\n",qm->SAstep,qm->SAsteps);
404 } /* write_gaussian_SH_input */
406 void write_gaussian_input(int step ,t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
416 out = fopen("input.com","w");
417 /* write the route */
419 if(qm->QMmethod>=eQMmethodRHF)
426 fprintf(out,"%s%3d\n",
427 "%nprocshare=",qm->nQMcpus);
428 fprintf(out,"%s%d\n",
430 /* use the modified links that include the LJ contribution at the QM level */
431 if(qm->bTS||qm->bOPT){
432 fprintf(out,"%s%s%s",
433 "%subst l701 ",qm->devel_dir,"/l701_LJ\n");
434 fprintf(out,"%s%s%s",
435 "%subst l301 ",qm->devel_dir,"/l301_LJ\n");
438 fprintf(out,"%s%s%s",
439 "%subst l701 ",qm->devel_dir,"/l701\n");
440 fprintf(out,"%s%s%s",
441 "%subst l301 ",qm->devel_dir,"/l301\n");
443 fprintf(out,"%s%s%s",
444 "%subst l9999 ",qm->devel_dir,"/l9999\n");
452 if(qm->QMmethod==eQMmethodB3LYPLAN){
454 "B3LYP/GEN Pseudo=Read");
458 eQMmethod_names[qm->QMmethod]);
460 if(qm->QMmethod>=eQMmethodRHF){
461 if(qm->QMmethod==eQMmethodCASSCF){
462 /* in case of cas, how many electrons and orbitals do we need?
464 fprintf(out,"(%d,%d)",
465 qm->CASelectrons,qm->CASorbitals);
468 eQMbasis_names[qm->QMbasis]);
471 if(QMMMrec->QMMMscheme==eQMMMschemenormal && mm->nrMMatoms){
475 if (step || qm->QMmethod==eQMmethodCASSCF){
476 /* fetch guess from checkpoint file, always for CASSCF */
477 fprintf(out,"%s"," guess=read");
479 fprintf(out,"\nNosymm units=bohr\n");
482 fprintf(out,"OPT=(Redundant,TS,noeigentest,ModRedundant) Punch=(Coord,Derivatives) ");
485 fprintf(out,"OPT=(Redundant,ModRedundant) Punch=(Coord,Derivatives) ");
488 fprintf(out,"FORCE Punch=(Derivatives) ");
490 fprintf(out,"iop(3/33=1)\n\n");
491 fprintf(out, "input-file generated by gromacs\n\n");
492 fprintf(out,"%2d%2d\n",qm->QMcharge,qm->multiplicity);
493 for (i=0;i<qm->nrQMatoms;i++){
495 fprintf(out,"%3d %10.7lf %10.7lf %10.7lf\n",
496 qm->atomicnumberQM[i],
497 qm->xQM[i][XX]/BOHR2NM,
498 qm->xQM[i][YY]/BOHR2NM,
499 qm->xQM[i][ZZ]/BOHR2NM);
501 fprintf(out,"%3d %10.7f %10.7f %10.7f\n",
502 qm->atomicnumberQM[i],
503 qm->xQM[i][XX]/BOHR2NM,
504 qm->xQM[i][YY]/BOHR2NM,
505 qm->xQM[i][ZZ]/BOHR2NM);
509 /* Pseudo Potential and ECP are included here if selected (MEthod suffix LAN) */
510 if(qm->QMmethod==eQMmethodB3LYPLAN){
512 for(i=0;i<qm->nrQMatoms;i++){
513 if(qm->atomicnumberQM[i]<21){
514 fprintf(out,"%d ",i+1);
517 fprintf(out,"\n%s\n****\n",eQMbasis_names[qm->QMbasis]);
519 for(i=0;i<qm->nrQMatoms;i++){
520 if(qm->atomicnumberQM[i]>21){
521 fprintf(out,"%d ",i+1);
524 fprintf(out,"\n%s\n****\n\n","lanl2dz");
526 for(i=0;i<qm->nrQMatoms;i++){
527 if(qm->atomicnumberQM[i]>21){
528 fprintf(out,"%d ",i+1);
531 fprintf(out,"\n%s\n","lanl2dz");
536 /* MM point charge data */
537 if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
538 fprintf(stderr,"nr mm atoms in gaussian.c = %d\n",mm->nrMMatoms);
540 if(qm->bTS||qm->bOPT){
541 /* freeze the frontier QM atoms and Link atoms. This is
542 * important only if a full QM subsystem optimization is done
543 * with a frozen MM environmeent. For dynamics, or gromacs's own
544 * optimization routines this is not important.
546 for(i=0;i<qm->nrQMatoms;i++){
547 if(qm->frontatoms[i]){
548 fprintf(out,"%d F\n",i+1); /* counting from 1 */
551 /* MM point charges include LJ parameters in case of QM optimization
553 for(i=0;i<mm->nrMMatoms;i++){
555 fprintf(out,"%10.7lf %10.7lf %10.7lf %8.4lf 0.0 %10.7lf %10.7lf\n",
556 mm->xMM[i][XX]/BOHR2NM,
557 mm->xMM[i][YY]/BOHR2NM,
558 mm->xMM[i][ZZ]/BOHR2NM,
560 mm->c6[i],mm->c12[i]);
562 fprintf(out,"%10.7f %10.7f %10.7f %8.4f 0.0 %10.7f %10.7f\n",
563 mm->xMM[i][XX]/BOHR2NM,
564 mm->xMM[i][YY]/BOHR2NM,
565 mm->xMM[i][ZZ]/BOHR2NM,
567 mm->c6[i],mm->c12[i]);
573 for(i=0;i<mm->nrMMatoms;i++){
575 fprintf(out,"%10.7lf %10.7lf %10.7lf %8.4lf\n",
576 mm->xMM[i][XX]/BOHR2NM,
577 mm->xMM[i][YY]/BOHR2NM,
578 mm->xMM[i][ZZ]/BOHR2NM,
581 fprintf(out,"%10.7f %10.7f %10.7f %8.4f\n",
582 mm->xMM[i][XX]/BOHR2NM,
583 mm->xMM[i][YY]/BOHR2NM,
584 mm->xMM[i][ZZ]/BOHR2NM,
595 } /* write_gaussian_input */
597 real read_gaussian_output(rvec QMgrad[],rvec MMgrad[],int step,
598 t_QMrec *qm, t_MMrec *mm)
609 in=fopen("fort.7","r");
613 /* in case of an optimization, the coordinates are printed in the
614 * fort.7 file first, followed by the energy, coordinates and (if
615 * required) the CI eigenvectors.
617 if(qm->bTS||qm->bOPT){
618 for(i=0;i<qm->nrQMatoms;i++){
619 if( NULL == fgets(buf,300,in))
621 gmx_fatal(FARGS,"Error reading Gaussian output - not enough atom lines?");
625 sscanf(buf,"%d %lf %lf %lf\n",
631 sscanf(buf,"%d %f %f %f\n",
638 qm->xQM[i][j]*=BOHR2NM;
642 /* the next line is the energy and in the case of CAS, the energy
643 * difference between the two states.
645 if(NULL == fgets(buf,300,in))
647 gmx_fatal(FARGS,"Error reading Gaussian output");
651 sscanf(buf,"%lf\n",&QMener);
653 sscanf(buf,"%f\n", &QMener);
655 /* next lines contain the gradients of the QM atoms */
656 for(i=0;i<qm->nrQMatoms;i++){
657 if(NULL == fgets(buf,300,in))
659 gmx_fatal(FARGS,"Error reading Gaussian output");
662 sscanf(buf,"%lf %lf %lf\n",
667 sscanf(buf,"%f %f %f\n",
673 /* the next lines are the gradients of the MM atoms */
674 if(qm->QMmethod>=eQMmethodRHF){
675 for(i=0;i<mm->nrMMatoms;i++){
676 if(NULL==fgets(buf,300,in))
678 gmx_fatal(FARGS,"Error reading Gaussian output");
681 sscanf(buf,"%lf %lf %lf\n",
686 sscanf(buf,"%f %f %f\n",
697 real read_gaussian_SH_output(rvec QMgrad[],rvec MMgrad[],int step,
698 gmx_bool swapped,t_QMrec *qm, t_MMrec *mm)
709 in=fopen("fort.7","r");
710 /* first line is the energy and in the case of CAS, the energy
711 * difference between the two states.
713 if(NULL == fgets(buf,300,in))
715 gmx_fatal(FARGS,"Error reading Gaussian output");
719 sscanf(buf,"%lf %lf\n",&QMener,&DeltaE);
721 sscanf(buf,"%f %f\n", &QMener,&DeltaE);
724 /* switch on/off the State Averaging */
726 if(DeltaE > qm->SAoff){
731 else if (DeltaE < qm->SAon || (qm->SAstep > 0)){
732 if (qm->SAstep < qm->SAsteps){
738 fprintf(stderr,"Gap = %5f,SA = %3d\n",DeltaE,(qm->SAstep>0));
739 /* next lines contain the gradients of the QM atoms */
740 for(i=0;i<qm->nrQMatoms;i++){
741 if(NULL==fgets(buf,300,in))
743 gmx_fatal(FARGS,"Error reading Gaussian output");
747 sscanf(buf,"%lf %lf %lf\n",
752 sscanf(buf,"%f %f %f\n",
758 /* the next lines, are the gradients of the MM atoms */
760 for(i=0;i<mm->nrMMatoms;i++){
761 if(NULL==fgets(buf,300,in))
763 gmx_fatal(FARGS,"Error reading Gaussian output");
766 sscanf(buf,"%lf %lf %lf\n",
771 sscanf(buf,"%f %f %f\n",
778 /* the next line contains the two CI eigenvector elements */
779 if(NULL==fgets(buf,300,in))
781 gmx_fatal(FARGS,"Error reading Gaussian output");
784 sscanf(buf,"%d",&qm->CIdim);
785 snew(qm->CIvec1,qm->CIdim);
786 snew(qm->CIvec1old,qm->CIdim);
787 snew(qm->CIvec2,qm->CIdim);
788 snew(qm->CIvec2old,qm->CIdim);
790 /* before reading in the new current CI vectors, copy the current
791 * CI vector into the old one.
793 for(i=0;i<qm->CIdim;i++){
794 qm->CIvec1old[i] = qm->CIvec1[i];
795 qm->CIvec2old[i] = qm->CIvec2[i];
799 for(i=0;i<qm->CIdim;i++){
800 if(NULL==fgets(buf,300,in))
802 gmx_fatal(FARGS,"Error reading Gaussian output");
805 sscanf(buf,"%lf\n",&qm->CIvec1[i]);
807 sscanf(buf,"%f\n", &qm->CIvec1[i]);
811 for(i=0;i<qm->CIdim;i++){
812 if(NULL==fgets(buf,300,in))
814 gmx_fatal(FARGS,"Error reading Gaussian output");
817 sscanf(buf,"%lf\n",&qm->CIvec2[i]);
819 sscanf(buf,"%f\n", &qm->CIvec2[i]);
826 real inproduct(real *a, real *b, int n)
833 /* computes the inner product between two vectors (a.b), both of
834 * which have length n.
842 int hop(int step, t_QMrec *qm)
847 d11=0.0,d12=0.0,d21=0.0,d22=0.0;
849 /* calculates the inproduct between the current Ci vector and the
850 * previous CI vector. A diabatic hop will be made if d12 and d21
851 * are much bigger than d11 and d22. In that case hop returns true,
852 * otherwise it returns false.
854 if(step){ /* only go on if more than one step has been done */
855 d11 = inproduct(qm->CIvec1,qm->CIvec1old,qm->CIdim);
856 d12 = inproduct(qm->CIvec1,qm->CIvec2old,qm->CIdim);
857 d21 = inproduct(qm->CIvec2,qm->CIvec1old,qm->CIdim);
858 d22 = inproduct(qm->CIvec2,qm->CIvec2old,qm->CIdim);
860 fprintf(stderr,"-------------------\n");
861 fprintf(stderr,"d11 = %13.8f\n",d11);
862 fprintf(stderr,"d12 = %13.8f\n",d12);
863 fprintf(stderr,"d21 = %13.8f\n",d21);
864 fprintf(stderr,"d22 = %13.8f\n",d22);
865 fprintf(stderr,"-------------------\n");
867 if((fabs(d12)>0.5)&&(fabs(d21)>0.5))
873 void do_gaussian(int step,char *exe)
878 /* make the call to the gaussian binary through system()
879 * The location of the binary will be picked up from the
880 * environment using getenv().
882 if(step) /* hack to prevent long inputfiles */
883 sprintf(buf,"%s < %s > %s",
888 sprintf(buf,"%s < %s > %s",
892 fprintf(stderr,"Calling '%s'\n",buf);
894 printf("Warning-- No calls to system(3) supported on this platform.");
895 gmx_fatal(FARGS,"Call to '%s' failed\n",buf);
897 if ( system(buf) != 0 )
898 gmx_fatal(FARGS,"Call to '%s' failed\n",buf);
902 real call_gaussian(t_commrec *cr, t_forcerec *fr,
903 t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
905 /* normal gaussian jobs */
918 sprintf(exe,"%s/%s",qm->gauss_dir,qm->gauss_exe);
919 snew(QMgrad,qm->nrQMatoms);
920 snew(MMgrad,mm->nrMMatoms);
922 write_gaussian_input(step,fr,qm,mm);
923 do_gaussian(step,exe);
924 QMener = read_gaussian_output(QMgrad,MMgrad,step,qm,mm);
925 /* put the QMMM forces in the force array and to the fshift
927 for(i=0;i<qm->nrQMatoms;i++){
929 f[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
930 fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
933 for(i=0;i<mm->nrMMatoms;i++){
935 f[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
936 fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
939 QMener = QMener*HARTREE2KJ*AVOGADRO;
944 } /* call_gaussian */
946 real call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
947 rvec f[], rvec fshift[])
949 /* a gaussian call routine intended for doing diabatic surface
950 * "sliding". See the manual for the theoretical background of this
960 swapped=FALSE; /* handle for identifying the current PES */
962 swap=FALSE; /* the actual swap */
971 sprintf(exe,"%s/%s",qm->gauss_dir,qm->gauss_exe);
972 /* hack to do ground state simulations */
975 buf = getenv("STATE");
977 sscanf(buf,"%d",&state);
986 /* copy the QMMMrec pointer */
987 snew(QMgrad,qm->nrQMatoms);
988 snew(MMgrad,mm->nrMMatoms);
989 /* at step 0 there should be no SA */
992 /* temporray set to step + 1, since there is a chk start */
993 write_gaussian_SH_input(step,swapped,fr,qm,mm);
995 do_gaussian(step,exe);
996 QMener = read_gaussian_SH_output(QMgrad,MMgrad,step,swapped,qm,mm);
998 /* check for a surface hop. Only possible if we were already state
1003 swap = (step && hop(step,qm));
1006 else { /* already on the other surface, so check if we go back */
1007 swap = (step && hop(step,qm));
1008 swapped =!swap; /* so swapped shoud be false again */
1010 if (swap){/* change surface, so do another call */
1011 write_gaussian_SH_input(step,swapped,fr,qm,mm);
1012 do_gaussian(step,exe);
1013 QMener = read_gaussian_SH_output(QMgrad,MMgrad,step,swapped,qm,mm);
1016 /* add the QMMM forces to the gmx force array and fshift
1018 for(i=0;i<qm->nrQMatoms;i++){
1020 f[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
1021 fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
1024 for(i=0;i<mm->nrMMatoms;i++){
1026 f[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
1027 fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
1030 QMener = QMener*HARTREE2KJ*AVOGADRO;
1031 fprintf(stderr,"step %5d, SA = %5d, swap = %5d\n",
1032 step,(qm->SAstep>0),swapped);
1037 } /* call_gaussian_SH */
1039 /* end of gaussian sub routines */
1043 gmx_qmmm_gaussian_empty;