1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
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-2008, 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 * Groningen Machine for Chemical Simulation
62 #include "gmx_fatal.h"
66 /* ORCA interface routines */
68 void init_orca(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
73 /* ORCA settings on the system */
74 buf = getenv("BASENAME");
76 snew(qm->orca_basename,200);
77 sscanf(buf,"%s",qm->orca_basename);
80 gmx_fatal(FARGS,"no $BASENAME\n");
82 /* ORCA directory on the system */
84 buf = getenv("ORCA_PATH");
85 fprintf(stderr,"%s",buf);
88 snew(qm->orca_dir,200);
89 sscanf(buf,"%s",qm->orca_dir);
92 gmx_fatal(FARGS,"no $ORCA_PATH, check manual\n");
94 fprintf(stderr,"%s...\n",qm->orca_dir);
95 fprintf(stderr,"orca initialised...\n");
96 /* since we append the output to the BASENAME.out file,
97 we should delete an existent old out-file here. */
98 sprintf(buf,"%s.out",qm->orca_basename);
103 void write_orca_input(int step ,t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
110 *out, *pcFile, *addInputFile, *LJCoeff;
112 *buf,*orcaInput,*addInputFilename,*LJCoeffFilename,
113 *pcFilename,*exclInName,*exclOutName;
115 /* write the first part of the input-file */
117 sprintf(orcaInput,"%s.inp",qm->orca_basename);
118 out = fopen(orcaInput,"w");
119 snew(addInputFilename,200);
120 sprintf(addInputFilename,"%s.ORCAINFO",qm->orca_basename);
121 addInputFile = fopen(addInputFilename,"r");
122 fprintf(out, "#input-file generated by gromacs\n");
124 fprintf(out,"!QMMMOpt TightSCF\n");
125 fprintf(out,"%s\n","%geom TS_Search EF end");
128 fprintf(out,"!QMMMOpt TightSCF\n");
131 fprintf(out,"!EnGrad TightSCF\n");
133 /* here we include the insertion of the additional orca-input */
135 if (addInputFile!=NULL) {
136 while (!feof(addInputFile)) {
137 if (fgets(buf, 200, addInputFile) != NULL)
142 fprintf(stderr,"No information on the calculation given in <%s>\n",addInputFilename);
143 gmx_call("qm_orca.c");
145 fclose(addInputFile);
146 if(qm->bTS||qm->bOPT){
147 /* freeze the frontier QM atoms and Link atoms. This is
148 * important only if a full QM subsystem optimization is done
149 * with a frozen MM environmeent. For dynamics, or gromacs's own
150 * optimization routines this is not important.
152 /* ORCA reads the exclusions from LJCoeffFilename.Excl,
153 *so we have to rename the file
156 for(i=0;i<qm->nrQMatoms;i++){
157 if(qm->frontatoms[i]){
159 fprintf(out,"%s\n","%geom");
160 fprintf(out," Constraints \n");
163 fprintf(out," {C %d C}\n",i); /* counting from 0 */
166 if (didStart) fprintf(out," end\n end\n");
167 /* make a file with information on the C6 and C12 coefficients */
168 if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
169 snew(exclInName,200);
170 snew(exclOutName,200);
171 sprintf(exclInName,"QMMMexcl.dat");
172 sprintf(exclOutName,"%s.LJ.Excl",qm->orca_basename);
173 rename(exclInName,exclOutName);
174 snew(LJCoeffFilename,200);
175 sprintf(LJCoeffFilename,"%s.LJ",qm->orca_basename);
176 fprintf(out,"%s%s%s\n","%LJCOEFFICIENTS \"",LJCoeffFilename,"\"");
177 /* make a file with information on the C6 and C12 coefficients */
178 LJCoeff = fopen(LJCoeffFilename,"w");
179 fprintf(LJCoeff,"%d\n",qm->nrQMatoms);
180 for (i=0;i<qm->nrQMatoms;i++){
182 fprintf(LJCoeff,"%10.7lf %10.7lf\n",qm->c6[i],qm->c12[i]);
184 fprintf(LJCoeff,"%10.7f %10.7f\n",qm->c6[i],qm->c12[i]);
187 fprintf(LJCoeff,"%d\n",mm->nrMMatoms);
188 for (i=0;i<mm->nrMMatoms;i++){
190 fprintf(LJCoeff,"%10.7lf %10.7lf\n",mm->c6[i],mm->c12[i]);
192 fprintf(LJCoeff,"%10.7f %10.7f\n",mm->c6[i],mm->c12[i]);
198 /* write charge and multiplicity
200 fprintf(out,"*xyz %2d%2d\n",qm->QMcharge,qm->multiplicity);
201 /* write the QM coordinates
203 for (i=0;i<qm->nrQMatoms;i++){
205 if (qm->atomicnumberQM[i]==0)
208 atomNr = qm->atomicnumberQM[i];
210 fprintf(out,"%3d %10.7lf %10.7lf %10.7lf\n",
216 fprintf(out,"%3d %10.7f %10.7f %10.7f\n",
224 /* write the MM point charge data
226 if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
227 /* name of the point charge file */
228 snew(pcFilename,200);
229 sprintf(pcFilename,"%s.pc",qm->orca_basename);
230 fprintf(out,"%s%s%s\n","%pointcharges \"",pcFilename,"\"");
231 pcFile = fopen(pcFilename,"w");
232 fprintf(pcFile,"%d\n",mm->nrMMatoms);
233 for(i=0;i<mm->nrMMatoms;i++){
235 fprintf(pcFile,"%8.4lf %10.7lf %10.7lf %10.7lf\n",
241 fprintf(pcFile,"%8.4f %10.7f %10.7f %10.7f\n",
248 fprintf(pcFile,"\n");
254 } /* write_orca_input */
256 real read_orca_output(rvec QMgrad[],rvec MMgrad[],int step,t_forcerec *fr,
257 t_QMrec *qm, t_MMrec *mm)
262 buf[300], tmp[300], orca_xyzFilename[300], orca_pcgradFilename[300], orca_engradFilename[300];
266 *xyz, *pcgrad, *engrad;
271 /* in case of an optimization, the coordinates are printed in the
272 * xyz file, the energy and gradients for the QM part are stored in the engrad file
273 * and the gradients for the point charges are stored in the pc file.
276 /* we need the new xyz coordinates of the QM atoms only for separate QM-optimization
279 if(qm->bTS||qm->bOPT){
280 sprintf(orca_xyzFilename,"%s.xyz",qm->orca_basename);
281 xyz=fopen(orca_xyzFilename,"r");
282 if (fgets(buf,300,xyz) == NULL)
283 gmx_fatal(FARGS, "Unexpected end of ORCA output");
284 if (fgets(buf,300,xyz) == NULL)
285 gmx_fatal(FARGS, "Unexpected end of ORCA output");
286 for(i=0;i<qm->nrQMatoms;i++){
287 if (fgets(buf,300,xyz) == NULL)
288 gmx_fatal(FARGS, "Unexpected end of ORCA output");
290 sscanf(buf,"%s%lf%lf%lf\n",
296 sscanf(buf,"%d%f%f%f\n",
308 sprintf(orca_engradFilename,"%s.engrad",qm->orca_basename);
309 engrad=fopen(orca_engradFilename,"r");
310 /* we read the energy and the gradient for the qm-atoms from the engrad file
312 /* we can skip the first seven lines
315 if (fgets(buf,300,engrad) == NULL)
316 gmx_fatal(FARGS, "Unexpected end of ORCA output");
318 /* now comes the energy
320 if (fgets(buf,300,engrad) == NULL)
321 gmx_fatal(FARGS, "Unexpected end of ORCA output");
323 sscanf(buf,"%lf\n",&QMener);
325 sscanf(buf,"%f\n", &QMener);
327 /* we can skip the next three lines
330 if (fgets(buf,300,engrad) == NULL)
331 gmx_fatal(FARGS, "Unexpected end of ORCA output");
333 /* next lines contain the gradients of the QM atoms
334 * now comes the gradient, one value per line:
335 * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
338 for(i=0;i<3*qm->nrQMatoms;i++){
340 if (fgets(buf,300,engrad) == NULL)
341 gmx_fatal(FARGS, "Unexpected end of ORCA output");
344 sscanf(buf,"%lf\n", &QMgrad[k][XX]);
346 sscanf(buf,"%lf\n", &QMgrad[k][YY]);
348 sscanf(buf,"%lf\n", &QMgrad[k][ZZ]);
351 sscanf(buf,"%f\n", &QMgrad[k][XX]);
353 sscanf(buf,"%f\n", &QMgrad[k][YY]);
355 sscanf(buf,"%f\n", &QMgrad[k][ZZ]);
359 /* write the MM point charge data
361 if(QMMMrec->QMMMscheme!=eQMMMschemeoniom && mm->nrMMatoms){
362 sprintf(orca_pcgradFilename,"%s.pcgrad",qm->orca_basename);
363 pcgrad=fopen(orca_pcgradFilename,"r");
365 /* we read the gradient for the mm-atoms from the pcgrad file
367 /* we can skip the first line
369 if (fgets(buf,300,pcgrad) == NULL)
370 gmx_fatal(FARGS, "Unexpected end of ORCA output");
371 for(i=0;i<mm->nrMMatoms;i++){
372 if (fgets(buf,300,pcgrad) == NULL)
373 gmx_fatal(FARGS, "Unexpected end of ORCA output");
375 sscanf(buf,"%lf%lf%lf\n",
380 sscanf(buf,"%f%f%f\n",
391 void do_orca(int step,char *exe, char *orca_dir, char *basename)
394 /* make the call to the orca binary through system()
395 * The location of the binary is set through the
400 sprintf(buf,"%s/%s %s.inp >> %s.out",
405 fprintf(stderr,"Calling '%s'\n",buf);
406 if ( system(buf) != 0 )
407 gmx_fatal(FARGS,"Call to '%s' failed\n",buf);
410 real call_orca(t_commrec *cr, t_forcerec *fr,
411 t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
413 /* normal orca jobs */
426 sprintf(exe,"%s","orca");
427 snew(QMgrad,qm->nrQMatoms);
428 snew(MMgrad,mm->nrMMatoms);
430 write_orca_input(step,fr,qm,mm);
431 do_orca(step,exe,qm->orca_dir,qm->orca_basename);
432 QMener = read_orca_output(QMgrad,MMgrad,step,fr,qm,mm);
433 /* put the QMMM forces in the force array and to the fshift
435 for(i=0;i<qm->nrQMatoms;i++){
437 f[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
438 fshift[i][j] = HARTREE_BOHR2MD*QMgrad[i][j];
441 for(i=0;i<mm->nrMMatoms;i++){
443 f[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
444 fshift[i+qm->nrQMatoms][j] = HARTREE_BOHR2MD*MMgrad[i][j];
447 QMener = QMener*HARTREE2KJ*AVOGADRO;
453 /* end of orca sub routines */