1 /* SemiEmpirical.c */
2 /**********************************************************************************************************
3 Copyright (c) 2002-2013 Abdul-Rahman Allouche. All rights reserved
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
6 documentation files (the Gabedit), to deal in the Software without restriction, including without limitation
7 the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
9
10 The above copyright notice and this permission notice shall be included in all copies or substantial portions
11 of the Software.
12
13 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
14 TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
15 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
16 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
17 DEALINGS IN THE SOFTWARE.
18 ************************************************************************************************************/
19 #include "../../Config.h"
20 #include <stdlib.h>
21 #include <stdio.h>
22 #include <math.h>
23 #ifndef G_OS_WIN32
24 #include <unistd.h>
25 #endif
26
27 #include "../Common/Global.h"
28 #include "../Utils/AtomsProp.h"
29 #include "../Utils/Utils.h"
30 #include "../Utils/Constants.h"
31 #include "../Geometry/Fragments.h"
32 #include "../Geometry/DrawGeom.h"
33 #include "../SemiEmpirical/AtomSE.h"
34 #include "../SemiEmpirical/MoleculeSE.h"
35 #include "../SemiEmpirical/SemiEmpiricalModel.h"
36 #include "../SemiEmpirical/SemiEmpirical.h"
37
38 static void calculateGradientMopac(SemiEmpiricalModel* seModel);
39 static void calculateEnergyMopac(SemiEmpiricalModel* seModel);
40 static void calculateGradientFireFly(SemiEmpiricalModel* seModel);
41 static void calculateEnergyFireFly(SemiEmpiricalModel* seModel);
42 static void calculateGradientOpenBabel(SemiEmpiricalModel* seModel);
43 static void calculateEnergyOpenBabel(SemiEmpiricalModel* seModel);
44 static void calculateGradientGeneric(SemiEmpiricalModel* seModel);
45 static void calculateEnergyGeneric(SemiEmpiricalModel* seModel);
46
47
48 /****************************************************************/
getMultiplicityName(gint multiplicity,gchar * buffer)49 static void getMultiplicityName(gint multiplicity, gchar* buffer)
50 {
51 if(multiplicity==1) sprintf(buffer,"Singlet");
52 else if(multiplicity==2) sprintf(buffer,"Doublet");
53 else if(multiplicity==3) sprintf(buffer,"Triplet");
54 else if(multiplicity==4) sprintf(buffer,"Quartet");
55 else if(multiplicity==5) sprintf(buffer,"Quintet");
56 else if(multiplicity==6) sprintf(buffer,"Sextet");
57 else sprintf(buffer,"UNKNOWN");
58 }
59 /*****************************************************************************/
getEnergyMopac(gchar * fileNameOut,gdouble * energy)60 static gboolean getEnergyMopac(gchar* fileNameOut, gdouble* energy)
61 {
62 FILE* file = NULL;
63 gchar buffer[1024];
64 gchar* pdest = NULL;
65
66 file = FOpen(fileNameOut, "rb");
67 if(!file) return FALSE;
68 while(!feof(file))
69 {
70 if(!fgets(buffer,BSIZE,file))break;
71 pdest = strstr( buffer, " FINAL HEAT OF FORMATION");
72 if(pdest)
73 {
74 pdest = strstr( buffer, "=");
75 if(pdest)
76 {
77 if(sscanf(pdest+1,"%lf",energy)==1)
78 {
79 fclose(file);
80 return TRUE;
81 }
82 }
83 }
84 }
85 fclose(file);
86 return FALSE;
87 }
88 /**********************************************************************************************************************************/
getGradientMopacOut(gchar * fileNameOut,SemiEmpiricalModel * seModel)89 static gboolean getGradientMopacOut(gchar* fileNameOut, SemiEmpiricalModel *seModel)
90 {
91 FILE* file = NULL;
92 gchar buffer[1024];
93 gchar* pdest = NULL;
94 gboolean Ok = FALSE;
95 gdouble tmp;
96 gint i;
97 gint j;
98
99 file = FOpen(fileNameOut, "rb");
100 /* fprintf(stderr,"DEBUG : Je suis dans getGradientMopac\n");*/
101 if(!file) return FALSE;
102 while(!feof(file))
103 {
104 if(!fgets(buffer,BSIZE,file))break;
105 if(strstr(buffer, "FINAL POINT AND DERIVATIVES"))break;
106 if(strstr(buffer, "FINAL") && strstr(buffer, "POINT") && strstr(buffer,"AND") && strstr(buffer,"DERIVATIVES"))break;
107 }
108 if(!strstr(buffer,"FINAL")) rewind(file);// Old Mopac, before 2018
109 while(!feof(file))
110 {
111 if(!fgets(buffer,BSIZE,file))break;
112 pdest = strstr( buffer, "PARAMETER ATOM TYPE VALUE GRADIENT");
113 if(pdest)
114 {
115 /* fprintf(stderr,"DEBUG : %s\n",pdest);*/
116 for(i=0;i<seModel->molecule.nAtoms;i++)
117 for(j=0;j<3;j++)
118 {
119 if(!fgets(buffer,BSIZE,file))break;
120 pdest = strstr( buffer, "CARTESIAN");
121 if(pdest)
122 {
123 /* fprintf(stderr,"DEBUG : %s\n",pdest+12);*/
124 if(sscanf(pdest+12,"%lf %lf",&tmp,&seModel->molecule.gradient[j][i])!=2)
125 {
126 fclose(file);
127 return FALSE;
128 }
129 }
130 /* fprintf(stderr,"DEBUG : %f\n",seModel->molecule.gradient[j][i]);*/
131 }
132 Ok = TRUE;
133 break;
134 }
135 pdest = strstr( buffer, "Cartesian Gradients"); /* MOZYME Keyword */
136 if(pdest)
137 {
138 gchar td[100];
139 gint d;
140 if(!fgets(buffer,BSIZE,file))break; /*Atom X ....*/
141 if(!fgets(buffer,BSIZE,file))break; /* backspace */
142 for(i=0;i<seModel->molecule.nAtoms;i++)
143 {
144 if(!fgets(buffer,BSIZE,file)) /* 1 O 0.000 -4.566 0.027 */
145 {
146 fclose(file);
147 return FALSE;
148 }
149 if(sscanf(buffer,"%d %s %lf %lf %lf",&d, td,
150 &seModel->molecule.gradient[0][i],
151 &seModel->molecule.gradient[1][i],
152 &seModel->molecule.gradient[2][i]
153 )
154 !=5)
155 {
156 fclose(file);
157 return FALSE;
158 }
159 }
160 Ok = TRUE;
161 break;
162 }
163 }
164 fclose(file);
165 return Ok;
166 }
167 /*****************************************************************************/
getGradientMopacAux(gchar * fileNameOut,SemiEmpiricalModel * seModel)168 static gboolean getGradientMopacAux(gchar* fileNameOut, SemiEmpiricalModel *seModel)
169 {
170 FILE* file = NULL;
171 gboolean Ok = FALSE;
172 gchar* fileNamePrefix = get_suffix_name_file(fileNameOut);
173 gchar* fileNameAux = g_strdup_printf("%s.aux", fileNamePrefix);
174 gchar** strlist = NULL;
175 gint n;
176
177 /*
178 fprintf(stderr,"fileNameOut %s\n",fileNameOut);
179 fprintf(stderr,"fileNamePrefix %s\n",fileNamePrefix);
180 fprintf(stderr,"fileNameAux %s\n",fileNameAux);
181 */
182 file = FOpen(fileNameAux, "rb");
183 if(file)
184 {
185 strlist = get_one_block_from_aux_mopac_file(file, "GRADIENTS:KCAL", &n);
186 if(strlist && n== seModel->molecule.nAtoms*3)
187 {
188 gint i;
189 gint j;
190 gint k;
191 k=0;
192 for(i=0;i<seModel->molecule.nAtoms;i++)
193 for(j=0;j<3;j++)
194 {
195 /* fprintf(stderr,"strlist %d %s\n",k,strlist[k]);*/
196 seModel->molecule.gradient[j][i] = atof(strlist[k]);
197 k++;
198 }
199 Ok = TRUE;
200 }
201 }
202 fclose(file);
203 if(strlist) strlist = free_one_string_table(strlist, n);
204 if(fileNamePrefix) g_free(fileNamePrefix);
205 if(fileNameAux) g_free(fileNameAux);
206 return Ok;
207 }
208 /*****************************************************************************/
getGradientMopac(gchar * fileNameOut,SemiEmpiricalModel * seModel)209 static gboolean getGradientMopac(gchar* fileNameOut, SemiEmpiricalModel *seModel)
210 {
211 gboolean Ok = FALSE;
212 /* try to read gradients from mopac out file : more digits than those of aux */
213 Ok = getGradientMopacOut(fileNameOut,seModel);
214 if(Ok) return Ok;
215 Ok = getGradientMopacAux(fileNameOut,seModel);
216 return Ok;
217 }
218 /*****************************************************************************/
runOneMopac(SemiEmpiricalModel * seModel,gchar * keyWords)219 static gchar* runOneMopac(SemiEmpiricalModel* seModel, gchar* keyWords)
220 {
221 FILE* file = NULL;
222 FILE* fileSH = NULL;
223 gint j;
224 gchar* fileNameIn = NULL;
225 gchar* fileNameOut = NULL;
226 gchar* fileNameSH = NULL;
227 gchar multiplicityStr[100];
228 gchar buffer[1024];
229 MoleculeSE m = seModel->molecule;
230 #ifdef G_OS_WIN32
231 gchar c='%';
232 #endif
233
234 if(m.nAtoms<1) return fileNameOut;
235 #ifndef G_OS_WIN32
236 fileNameSH = g_strdup_printf("%s%sMopacOne.sh",seModel->workDir,G_DIR_SEPARATOR_S);
237 #else
238 fileNameSH = g_strdup_printf("%s%sMopacOne.bat",seModel->workDir,G_DIR_SEPARATOR_S);
239 #endif
240 fileSH = FOpen(fileNameSH, "w");
241 if(!fileSH) return FALSE;
242 #ifdef G_OS_WIN32
243 fprintf(fileSH,"@echo off\n");
244 fprintf(fileSH,"set PATH=%cPATH%c;\"%s\"\n",c,c,mopacDirectory);
245 #endif
246
247 getMultiplicityName(seModel->molecule.spinMultiplicity, multiplicityStr);
248
249 fileNameIn = g_strdup_printf("%s%sOne.mop",seModel->workDir,G_DIR_SEPARATOR_S);
250 file = FOpen(fileNameIn, "w");
251 if(!file)
252 {
253 if(fileNameIn) g_free(fileNameIn);
254 if(fileNameOut) g_free(fileNameOut);
255 if(fileNameSH) g_free(fileNameSH);
256 return FALSE;
257 }
258 fprintf(file,"* ===============================\n");
259 fprintf(file,"* Input file for Mopac\n");
260 fprintf(file,"* ===============================\n");
261 if(seModel->molecule.spinMultiplicity>1)
262 fprintf(file,"%s UHF CHARGE=%d AUX %s\n",keyWords,seModel->molecule.totalCharge,multiplicityStr);
263 else
264 fprintf(file,"%s CHARGE=%d AUX %s\n",keyWords,seModel->molecule.totalCharge,multiplicityStr);
265 fprintf(file,"\n");
266 fprintf(file,"Mopac file generated by Gabedit\n");
267
268 for(j=0;j<m.nAtoms;j++)
269 {
270 fprintf(file," %s %f %d %f %d %f %d\n",
271 m.atoms[j].prop.symbol,
272 m.atoms[j].coordinates[0],
273 1,
274 m.atoms[j].coordinates[1],
275 1,
276 m.atoms[j].coordinates[2],
277 1
278 );
279 }
280 fclose(file);
281 #ifndef G_OS_WIN32
282 fprintf(fileSH,"%s %s\n",NameCommandMopac,fileNameIn);
283 fclose(fileSH);
284 sprintf(buffer,"chmod u+x %s",fileNameSH);
285 {int ierr = system(buffer);}
286 { int ierr = system(fileNameSH);}
287 #else
288 fprintf(fileSH,"\"%s\" \"%s\"\n",NameCommandMopac,fileNameIn);
289 fclose(fileSH);
290 sprintf(buffer,"\"%s\"",fileNameSH);
291 {int ierr = system(buffer);}
292 #endif
293
294 unlink(fileNameIn);
295 unlink(fileNameSH);
296 if(fileNameIn) g_free(fileNameIn);
297 if(fileNameSH) g_free(fileNameSH);
298 fileNameOut = g_strdup_printf("%s%sOne.out",seModel->workDir,G_DIR_SEPARATOR_S);
299 return fileNameOut;
300 }
301 /**********************************************************************/
newMopacModel(gchar * method,gchar * dirName,SemiEmpiricalModelConstraints constraints)302 static SemiEmpiricalModel newMopacModel(gchar* method, gchar* dirName, SemiEmpiricalModelConstraints constraints)
303 {
304 SemiEmpiricalModel seModel = newSemiEmpiricalModel(method, dirName, constraints);
305
306 seModel.klass->calculateGradient = calculateGradientMopac;
307 seModel.klass->calculateEnergy = calculateEnergyMopac;
308
309 return seModel;
310 }
311 /**********************************************************************/
calculateGradientMopac(SemiEmpiricalModel * seModel)312 static void calculateGradientMopac(SemiEmpiricalModel* seModel)
313 {
314 gint i;
315 gint j;
316 MoleculeSE m = seModel->molecule;
317 gchar* keyWords = NULL;
318 gchar* fileOut = NULL;
319 if(!seModel) return;
320 if(seModel->molecule.nAtoms<1) return;
321 if(!seModel->method) return;
322 keyWords = g_strdup_printf("%s 1SCF GRAD",seModel->method);
323 fileOut = runOneMopac(seModel, keyWords);
324
325 if(fileOut)
326 {
327 for(j=0;j<3;j++)
328 for( i=0; i<m.nAtoms;i++)
329 m.gradient[j][i] = 0.0;
330 if(!getGradientMopac(fileOut, seModel))
331 {
332 StopCalcul=TRUE;
333 set_text_to_draw(_("Problem : I cannot caculate the Gradient... "));
334 set_statubar_operation_str(_("Calculation Stopped "));
335 drawGeom();
336 gtk_widget_set_sensitive(StopButton, FALSE);
337 Waiting(1);
338 return;
339 }
340 getEnergyMopac(fileOut, &seModel->molecule.energy);
341 g_free(fileOut);
342 }
343
344 }
345 /**********************************************************************/
calculateEnergyMopac(SemiEmpiricalModel * seModel)346 static void calculateEnergyMopac(SemiEmpiricalModel* seModel)
347 {
348 gchar* keyWords = NULL;
349 gchar* fileOut = NULL;
350 if(!seModel) return;
351 if(seModel->molecule.nAtoms<1) return;
352 if(!seModel->method) return;
353 keyWords = g_strdup_printf("%s 1SCF",seModel->method);
354 fileOut = runOneMopac(seModel, keyWords);
355 if(fileOut)
356 {
357 getEnergyMopac(fileOut, &seModel->molecule.energy);
358 g_free(fileOut);
359 }
360
361 }
362 /**********************************************************************/
createMopacModel(GeomDef * geom,gint Natoms,gint charge,gint spin,gchar * method,gchar * dirName,SemiEmpiricalModelConstraints constraints)363 SemiEmpiricalModel createMopacModel (GeomDef* geom,gint Natoms,gint charge, gint spin, gchar* method, gchar* dirName, SemiEmpiricalModelConstraints constraints)
364 {
365 SemiEmpiricalModel seModel = newMopacModel(method, dirName, constraints);
366
367 seModel.molecule = createMoleculeSE(geom,Natoms, charge, spin,TRUE);
368 setRattleConstraintsParameters(&seModel);
369
370 return seModel;
371 }
372 /**********************************************************************/
getEnergyFireFly(gchar * fileNameOut,gdouble * energy)373 static gboolean getEnergyFireFly(gchar* fileNameOut, gdouble* energy)
374 {
375 FILE* file = NULL;
376 gchar buffer[1024];
377 gchar* pdest = NULL;
378 gboolean OK = FALSE;
379
380 file = FOpen(fileNameOut, "rb");
381 if(!file) return FALSE;
382 while(!feof(file))
383 {
384 if(!fgets(buffer,BSIZE,file))break;
385 pdest = strstr( buffer, "HEAT OF FORMATION IS");
386 if(pdest)
387 {
388 pdest = strstr( buffer, "S");
389 if(pdest)
390 {
391 if(sscanf(pdest+1,"%lf",energy)==1)
392 {
393 OK = TRUE;
394 /* break;*/
395 }
396 }
397 }
398 }
399 fclose(file);
400 return OK;
401 }
402 /*****************************************************************************/
getGradientFireFly(gchar * fileNameOut,SemiEmpiricalModel * seModel)403 static gboolean getGradientFireFly(gchar* fileNameOut, SemiEmpiricalModel *seModel)
404 {
405 FILE* file = NULL;
406 gchar buffer[1024];
407 gchar stmp[1024];
408 gchar* pdest = NULL;
409 gboolean Ok = FALSE;
410 gint itmp;
411 gint i;
412 gint j;
413
414 file = FOpen(fileNameOut, "rb");
415 if(!file) return FALSE;
416 while(!feof(file))
417 {
418 if(!fgets(buffer,BSIZE,file))break;
419 pdest = strstr( buffer, "ATOM E'X E'Y E'Z");
420 if(pdest)
421 {
422 for(i=0;i<seModel->molecule.nAtoms;i++)
423 {
424 if(!fgets(buffer,BSIZE,file))break;
425 if(sscanf(buffer,"%d %s %lf %lf %lf",&itmp, stmp,
426 &seModel->molecule.gradient[0][i],
427 &seModel->molecule.gradient[1][i],
428 &seModel->molecule.gradient[2][i]
429 )!=5)
430 {
431 fclose(file);
432 return FALSE;
433 }
434 for(j=0;j<3;j++) seModel->molecule.gradient[j][i] *= 627.50944796/BOHR_TO_ANG;
435 }
436 Ok = TRUE;
437 break;
438 }
439 }
440 fclose(file);
441 return Ok;
442 }
443 /*****************************************************************************/
runOneFireFly(SemiEmpiricalModel * seModel,gchar * keyWords)444 static gchar* runOneFireFly(SemiEmpiricalModel* seModel, gchar* keyWords)
445 {
446 FILE* file = NULL;
447 FILE* fileSH = NULL;
448 gint j;
449 gchar* fileNameIn = NULL;
450 gchar* fileNameOut = NULL;
451 gchar* fileNameSH = NULL;
452 gchar multiplicityStr[100];
453 gchar buffer[1024];
454 MoleculeSE m = seModel->molecule;
455 gchar* fileNamePrefix = NULL;
456 #ifdef G_OS_WIN32
457 gchar c='%';
458 #endif
459
460 if(m.nAtoms<1) return fileNameOut;
461 #ifndef G_OS_WIN32
462 fileNameSH = g_strdup_printf("%s%sFireFlyOne.sh",seModel->workDir,G_DIR_SEPARATOR_S);
463 #else
464 fileNameSH = g_strdup_printf("%s%sFireFlyOne.bat",seModel->workDir,G_DIR_SEPARATOR_S);
465 #endif
466 fileSH = FOpen(fileNameSH, "w");
467 if(!fileSH) return FALSE;
468 #ifdef G_OS_WIN32
469 fprintf(fileSH,"@echo off\n");
470 fprintf(fileSH,"set PATH=%cPATH%c;\"%s\"\n",c,c,fireflyDirectory);
471 #endif
472
473 getMultiplicityName(seModel->molecule.spinMultiplicity, multiplicityStr);
474
475 fileNameIn = g_strdup_printf("%s%sOne.inp",seModel->workDir,G_DIR_SEPARATOR_S);
476 fileNameOut = g_strdup_printf("%s%sOne.out",seModel->workDir,G_DIR_SEPARATOR_S);
477
478
479 file = FOpen(fileNameIn, "w");
480 if(!file)
481 {
482 if(fileNameIn) g_free(fileNameIn);
483 if(fileNameOut) g_free(fileNameOut);
484 if(fileNameSH) g_free(fileNameSH);
485 return FALSE;
486 }
487 fprintf(file,"! ======================================================\n");
488 fprintf(file,"! Input file for FireFly\n");
489 fprintf(file,"! ======================================================\n");
490 if(strstr(keyWords,"RUNTYP"))
491 {
492 sscanf(strstr(keyWords,"RUNTYP"),"%s",buffer);
493 fprintf(file," $CONTRL %s $END\n",buffer);
494 }
495 if(strstr(keyWords,"SCFTYP"))
496 {
497 sscanf(strstr(keyWords,"SCFTYP"),"%s",buffer);
498 fprintf(file," $CONTRL %s $END\n",buffer);
499 }
500 else
501 {
502 if(seModel->molecule.spinMultiplicity==1)
503 fprintf(file," $CONTRL SCFTYP=RHF $END\n");
504 else
505 fprintf(file," $CONTRL SCFTYP=UHF $END\n");
506 }
507
508 fprintf(file," $CONTRL ICHARG=%d MULT=%d $END\n",seModel->molecule.totalCharge,seModel->molecule.spinMultiplicity);
509 if(strstr(keyWords,"GBASIS"))
510 {
511 sscanf(strstr(keyWords,"GBASIS"),"%s",buffer);
512 fprintf(file," $BASIS %s $END\n",buffer);
513 }
514 fprintf(file," $DATA\n");
515 fprintf(file,"Molecule specification\n");
516 fprintf(file,"C1\n");
517 for(j=0;j<m.nAtoms;j++)
518 {
519 gchar* symbol = m.atoms[j].prop.symbol;
520 SAtomsProp prop = prop_atom_get(symbol);
521 fprintf(file,"%s %f %f %f %f\n",
522 symbol,
523 (gdouble)prop.atomicNumber,
524 m.atoms[j].coordinates[0],
525 m.atoms[j].coordinates[1],
526 m.atoms[j].coordinates[2]
527 );
528 }
529 fprintf(file," $END\n");
530 fclose(file);
531 fileNamePrefix = g_strdup_printf("%s%sWorkFF",seModel->workDir,G_DIR_SEPARATOR_S);
532 #ifndef G_OS_WIN32
533 if(!strcmp(NameCommandFireFly,"pcgamess") || !strcmp(NameCommandFireFly,"nohup pcgamess")
534 || !strcmp(NameCommandFireFly,"firefly") || !strcmp(NameCommandFireFly,"nohup firefly"))
535 {
536 fprintf(fileSH,"mkdir %stmp\n",fileNamePrefix);
537 fprintf(fileSH,"cd %stmp\n",fileNamePrefix);
538 fprintf(fileSH,"cp %s input\n",fileNameIn);
539 fprintf(fileSH,"%s -p -o %s\n",NameCommandFireFly,fileNameOut);
540 fprintf(fileSH,"cd ..\n");
541 fprintf(fileSH,"rm PUNCH\n");
542 fprintf(fileSH,"/bin/rm -r %stmp\n",fileNamePrefix);
543 }
544 else
545 fprintf(fileSH,"%s %s",NameCommandFireFly,fileNameIn);
546 #else
547 if(!strcmp(NameCommandFireFly,"pcgamess") ||
548 !strcmp(NameCommandFireFly,"firefly") )
549 {
550 fprintf(fileSH,"mkdir \"%stmp\"\n",fileNamePrefix);
551 addUnitDisk(fileSH, fileNamePrefix);
552 fprintf(fileSH,"cd \"%stmp\"\n",fileNamePrefix);
553 fprintf(fileSH,"copy \"%s\" input\n",fileNameIn);
554 fprintf(fileSH,"%s -p -o \"%s\"\n",NameCommandFireFly,fileNameOut);
555 fprintf(fileSH,"cd ..\n");
556 fprintf(fileSH,"del PUNCH 2>nul\n");
557 fprintf(fileSH,"del /Q \"%stmp\"\n",fileNamePrefix);
558 fprintf(fileSH,"rmdir \"%stmp\"\n",fileNamePrefix);
559 }
560 else
561 fprintf(fileSH,"%s %s",NameCommandFireFly,fileNameIn);
562 #endif
563 fclose(fileSH);
564 #ifndef G_OS_WIN32
565 sprintf(buffer,"chmod u+x %s",fileNameSH);
566 {int ierr = system(buffer);}
567 { int ierr = system(fileNameSH);}
568 #else
569 sprintf(buffer,"\"%s\"",fileNameSH);
570 {int ierr = system(buffer);}
571 #endif
572 unlink(fileNameIn);
573 unlink(fileNameSH);
574 if(fileNamePrefix) g_free(fileNamePrefix);
575 if(fileNameIn) g_free(fileNameIn);
576 if(fileNameSH) g_free(fileNameSH);
577 return fileNameOut;
578 }
579 /**********************************************************************/
newFireFlyModel(gchar * method,gchar * dirName,SemiEmpiricalModelConstraints constraints)580 static SemiEmpiricalModel newFireFlyModel(gchar* method, gchar* dirName, SemiEmpiricalModelConstraints constraints)
581 {
582 SemiEmpiricalModel seModel = newSemiEmpiricalModel(method, dirName, constraints);
583
584 seModel.klass->calculateGradient = calculateGradientFireFly;
585 seModel.klass->calculateEnergy = calculateEnergyFireFly;
586
587 return seModel;
588 }
589 /**********************************************************************/
calculateGradientFireFly(SemiEmpiricalModel * seModel)590 static void calculateGradientFireFly(SemiEmpiricalModel* seModel)
591 {
592 gint i;
593 gint j;
594 MoleculeSE m = seModel->molecule;
595 gchar* keyWords = NULL;
596 gchar* fileOut = NULL;
597 if(!seModel) return;
598 if(seModel->molecule.nAtoms<1) return;
599 if(!seModel->method) return;
600 keyWords = g_strdup_printf("RUNTYP=GRADIENT GBASIS=%s",seModel->method);
601 fileOut = runOneFireFly(seModel, keyWords);
602
603 if(fileOut)
604 {
605 for(j=0;j<3;j++)
606 for( i=0; i<m.nAtoms;i++)
607 m.gradient[j][i] = 0.0;
608 if(!getGradientFireFly(fileOut, seModel))
609 {
610 #ifdef G_OS_WIN32
611 gchar* comm = g_strdup_printf("type %s",fileOut);
612 #else
613 gchar* comm = g_strdup_printf("cat %s",fileOut);
614 #endif
615 StopCalcul=TRUE;
616 set_text_to_draw(_("Problem : I cannot caculate the Gradient... "));
617 set_statubar_operation_str(_("Calculation Stopped "));
618 drawGeom();
619 gtk_widget_set_sensitive(StopButton, FALSE);
620 Waiting(1);
621 {int ierr = system(comm);}
622 g_free(fileOut);
623 g_free(comm);
624 return;
625 }
626 getEnergyFireFly(fileOut, &seModel->molecule.energy);
627 g_free(fileOut);
628 }
629
630 }
631 /**********************************************************************/
calculateEnergyFireFly(SemiEmpiricalModel * seModel)632 static void calculateEnergyFireFly(SemiEmpiricalModel* seModel)
633 {
634 gchar* keyWords = NULL;
635 gchar* fileOut = NULL;
636 if(!seModel) return;
637 if(seModel->molecule.nAtoms<1) return;
638 if(!seModel->method) return;
639 keyWords = g_strdup_printf("RUNTYP=Energy GBASIS=%s",seModel->method);
640 fileOut = runOneFireFly(seModel, keyWords);
641 if(fileOut)
642 {
643 getEnergyFireFly(fileOut, &seModel->molecule.energy);
644 g_free(fileOut);
645 }
646
647 }
648 /**********************************************************************/
createFireFlyModel(GeomDef * geom,gint Natoms,gint charge,gint spin,gchar * method,gchar * dirName,SemiEmpiricalModelConstraints constraints)649 SemiEmpiricalModel createFireFlyModel (GeomDef* geom,gint Natoms,gint charge, gint spin, gchar* method, gchar* dirName, SemiEmpiricalModelConstraints constraints)
650 {
651 SemiEmpiricalModel seModel = newFireFlyModel(method,dirName, constraints);
652
653 seModel.molecule = createMoleculeSE(geom,Natoms, charge, spin,TRUE);
654 setRattleConstraintsParameters(&seModel);
655
656 return seModel;
657 }
658 /**********************************************************************/
getEnergyOpenBabel(gchar * fileNameOut,gdouble * energy)659 static gboolean getEnergyOpenBabel(gchar* fileNameOut, gdouble* energy)
660 {
661 FILE* file = NULL;
662 char buffer[1024];
663 char* pdest = NULL;
664 //char* energyTag = "TOTAL ENERGY =";
665 char* energyTag = "FINAL ENERGY:";
666
667 file = fopen(fileNameOut, "r");
668 if(!file) return FALSE;
669 while(!feof(file))
670 {
671 if(!fgets(buffer,BSIZE,file))break;
672 pdest = strstr( buffer, energyTag);
673 if(pdest &&sscanf(pdest+strlen(energyTag)+1,"%lf",energy)==1)
674 {
675 fclose(file);
676 if(strstr(pdest,"kJ")) *energy /= KCALTOKJ;
677 return TRUE;
678 }
679 }
680 fclose(file);
681 return FALSE;
682 }
683 /*****************************************************************************/
getGradientOpenBabel(gchar * fileNameOut,SemiEmpiricalModel * seModel)684 static gboolean getGradientOpenBabel(gchar* fileNameOut, SemiEmpiricalModel *seModel)
685 {
686 FILE* file = NULL;
687 char buffer[1024];
688 char* pdest = NULL;
689 //char* energyTag = "TOTAL ENERGY =";
690 char* energyTag = "FINAL ENERGY:";
691 char* gradTag = "Gradients:";
692 gboolean kj = FALSE;
693 double conv = 1.0;
694 int i;
695 MoleculeSE* mol = &seModel->molecule;
696
697 file = fopen(fileNameOut, "r");
698 if(!file) return FALSE;
699 while(!feof(file))
700 {
701 if(!fgets(buffer,BSIZE,file))break;
702 pdest = strstr( buffer, energyTag);
703 if(pdest &&sscanf(pdest+strlen(energyTag)+1,"%lf",&mol->energy)==1)
704 {
705 if(strstr(pdest,"kJ")) { kj = TRUE;}
706 break;
707 }
708 }
709 while(!feof(file))
710 {
711 if(!fgets(buffer,BSIZE,file))break;
712 if(strstr(buffer, gradTag))
713 {
714 for(i=0;i<mol->nAtoms;i++)
715 {
716 if(!fgets(buffer,BSIZE,file))break;
717 //printf("%s\n",buffer);
718 if(sscanf(buffer,"%lf %lf %lf",
719 &mol->gradient[0][i],
720 &mol->gradient[1][i],
721 &mol->gradient[2][i]
722 )!=3) break;
723 }
724 break;
725 }
726 }
727 if(kj) conv /= KCALTOKJ;
728 mol->energy *= conv;
729 for(i=0;i<mol->nAtoms;i++)
730 {
731 mol->gradient[0][i] *= conv;
732 mol->gradient[1][i] *= conv;
733 mol->gradient[2][i] *= conv;
734 }
735
736 fclose(file);
737 return FALSE;
738 }
739 /*****************************************************************************/
runOneOpenBabel(SemiEmpiricalModel * seModel,gchar * keyWords)740 static gchar* runOneOpenBabel(SemiEmpiricalModel* seModel, gchar* keyWords)
741 {
742 FILE* fileSH = NULL;
743 char* fileNameIn = NULL;
744 char* fileNameOut = NULL;
745 char* fileNameSH = NULL;
746 char buffer[1024];
747 MoleculeSE* mol = &seModel->molecule;
748 char* NameCommandOpenBabel = "obgradient";
749 int rank = 0;
750 #ifdef ENABLE_MPI
751 MPI_Comm_rank( MPI_COMM_WORLD,&rank);
752 #endif
753 #ifdef G_OS_WIN32
754 char c='%';
755 #endif
756
757 if(mol->nAtoms<1) return fileNameOut;
758 #ifndef G_OS_WIN32
759 fileNameSH = g_strdup_printf("%s%sOpenBabelOne%d.sh",seModel->workDir,G_DIR_SEPARATOR_S,rank);
760 #else
761 fileNameSH = g_strdup_printf("%s%sOpenBabelOne%d.bat",seModel->workDir,G_DIR_SEPARATOR_S,rank);
762 #endif
763 fileSH = fopen(fileNameSH, "w");
764 if(!fileSH) return FALSE;
765 #ifdef G_OS_WIN32
766 fprintf(fileSH,"@echo off\n");
767 #endif
768
769 fileNameIn = g_strdup_printf("%s%sOne%d.hin",seModel->workDir,G_DIR_SEPARATOR_S,rank);
770 fileNameOut = g_strdup_printf("%s%sOne%d.out",seModel->workDir,G_DIR_SEPARATOR_S,rank);
771
772 if(!saveMoleculeSEHIN(mol,fileNameIn))
773 {
774 if(fileNameIn) free(fileNameIn);
775 if(fileNameOut) free(fileNameOut);
776 if(fileNameSH) free(fileNameSH);
777 return FALSE;
778 }
779 #ifndef G_OS_WIN32
780 if(keyWords)
781 {
782 fprintf(fileSH,"#!/bin/bash\n");
783 fprintf(fileSH,"export PATH=$PATH:%s\n",openbabelDirectory);
784 fprintf(fileSH,"export BABEL_DATADIR=%s\n",openbabelDirectory);
785 fprintf(fileSH,"%s %s %s > %s 2>/dev/null\n",NameCommandOpenBabel,keyWords,fileNameIn,fileNameOut);
786 fprintf(fileSH,"exit\n");
787 }
788 else
789 fprintf(fileSH,"%s %s >%s 2>/dev/null",NameCommandOpenBabel,fileNameIn, fileNameOut);
790 #else
791 {
792 if(strstr(openbabelDirectory,"\""))
793 {
794 fprintf(fileSH,"set PATH=%s;%cPATH%c\n",openbabelDirectory,'%','%');
795 fprintf(fileSH,"set BABEL_DATADIR=%s\n",openbabelDirectory);
796 }
797 else
798 {
799 fprintf(fileSH,"set PATH=\"%s\";%cPATH%c\n",openbabelDirectory,'%','%');
800 fprintf(fileSH,"set BABEL_DATADIR=%s\n",openbabelDirectory);
801 }
802
803 fprintf(fileSH,"%s %s %s > %s\n",NameCommandOpenBabel,keyWords,fileNameIn,fileNameOut);
804 fprintf(fileSH,"exit\n");
805 }
806 #endif
807 fclose(fileSH);
808 #ifndef G_OS_WIN32
809
810 /*
811 sprintf(buffer,"cat %s",fileNameSH);
812 system(buffer);
813 sprintf(buffer,"cat %s",fileNameIn);
814 system(buffer);
815 */
816
817 sprintf(buffer,"chmod u+x %s",fileNameSH);
818 system(buffer);
819 system(fileNameSH);
820
821 /*
822 sprintf(buffer,"cat %s",fileNameOut);
823 system(buffer);
824 */
825 #else
826 sprintf(buffer,"\"%s\"",fileNameSH);
827 system(buffer);
828 #endif
829 unlink(fileNameIn);
830 unlink(fileNameSH);
831 if(fileNameIn) free(fileNameIn);
832 if(fileNameSH) free(fileNameSH);
833 return fileNameOut;
834 }
835 /**********************************************************************/
newOpenBabelModel(gchar * method,gchar * dirName,SemiEmpiricalModelConstraints constraints)836 static SemiEmpiricalModel newOpenBabelModel(gchar* method, gchar* dirName, SemiEmpiricalModelConstraints constraints)
837 {
838 SemiEmpiricalModel seModel = newSemiEmpiricalModel(method, dirName, constraints);
839
840 seModel.klass->calculateGradient = calculateGradientOpenBabel;
841 seModel.klass->calculateEnergy = calculateEnergyOpenBabel;
842
843 return seModel;
844 }
845 /**********************************************************************/
calculateGradientOpenBabel(SemiEmpiricalModel * seModel)846 static void calculateGradientOpenBabel(SemiEmpiricalModel* seModel)
847 {
848 char* keyWords = NULL;
849 char* fileOut = NULL;
850 if(!seModel) return;
851 if(seModel->molecule.nAtoms<1) return;
852 if(!seModel->method) return;
853 keyWords = g_strdup_printf("-ff %s",seModel->method);
854 fileOut = runOneOpenBabel(seModel, keyWords);
855 if(fileOut)
856 {
857 getGradientOpenBabel(fileOut, seModel);
858 //getDipoleOpenBabel(fileOut, &seModel->molecule, seModel->molecule.dipole);
859 computeMoleculeSEDipole(&seModel->molecule);
860 free(fileOut);
861 }
862
863 }
864 /**********************************************************************/
calculateEnergyOpenBabel(SemiEmpiricalModel * seModel)865 static void calculateEnergyOpenBabel(SemiEmpiricalModel* seModel)
866 {
867 char* keyWords = NULL;
868 char* fileOut = NULL;
869 if(!seModel) return;
870 if(seModel->molecule.nAtoms<1) return;
871 if(!seModel->method) return;
872 keyWords = g_strdup_printf("%s ",seModel->method);
873 fileOut = runOneOpenBabel(seModel, keyWords);
874 if(fileOut)
875 {
876 getEnergyOpenBabel(fileOut, &seModel->molecule.energy);
877 //getDipoleOpenBabel(fileOut, &seModel->molecule, seModel->molecule.dipole);
878 computeMoleculeSEDipole(&seModel->molecule);
879 free(fileOut);
880 }
881
882
883 }
884 /**********************************************************************/
createOpenBabelModel(GeomDef * geom,gint Natoms,gint charge,gint spin,gchar * method,gchar * dirName,SemiEmpiricalModelConstraints constraints)885 SemiEmpiricalModel createOpenBabelModel (GeomDef* geom,gint Natoms,gint charge, gint spin, gchar* method, gchar* dirName, SemiEmpiricalModelConstraints constraints)
886 {
887 SemiEmpiricalModel seModel = newOpenBabelModel(method,dirName, constraints);
888
889 seModel.molecule = createMoleculeSE(geom,Natoms, charge, spin,TRUE);
890 setRattleConstraintsParameters(&seModel);
891
892 return seModel;
893 }
894 /**********************************************************************/
getDipoleGeneric(char * fileNameOut,double * dipole)895 static gboolean getDipoleGeneric(char* fileNameOut, double* dipole)
896 {
897 FILE* file = NULL;
898 char buffer[1024];
899 int i;
900 file = fopen(fileNameOut, "r");
901 if(!file) return FALSE;
902 if(!fgets(buffer,BSIZE,file)) { fclose(file); return FALSE;}/* first line for energy in Hartree*/
903 if(!fgets(buffer,BSIZE,file)) { fclose(file); return FALSE;}/* dipole in au */
904 for(i=0;i<strlen(buffer);i++) if(buffer[i]=='D' || buffer[i]=='d') buffer[i] ='E';
905 if(sscanf(buffer,"%lf %lf %lf",&dipole[0],&dipole[1],&dipole[2])==3)
906 {
907 for(i=0;i<3;i++) dipole[i] *= AUTODEB;
908 fclose(file);
909 return TRUE;
910 }
911 fclose(file);
912 return FALSE;
913 }
914 /*****************************************************************************/
getEnergyGeneric(char * fileNameOut,double * energy)915 static gboolean getEnergyGeneric(char* fileNameOut, double* energy)
916 {
917 FILE* file = NULL;
918 char buffer[1024];
919 int i;
920 file = fopen(fileNameOut, "r");
921 if(!file) return FALSE;
922 if(!fgets(buffer,BSIZE,file)) { fclose(file); return FALSE;}/* first line for energy in Hartree*/
923
924 for(i=0;i<strlen(buffer);i++) if(buffer[i]=='D' || buffer[i]=='d') buffer[i] ='E';
925 if(sscanf(buffer,"%lf",energy)==1)
926 {
927 fclose(file);
928 *energy *= AUTOKCAL;
929 return TRUE;
930 }
931 fclose(file);
932 return FALSE;
933 }
934 /*****************************************************************************/
getGradientGeneric(char * fileNameOut,SemiEmpiricalModel * seModel)935 gboolean getGradientGeneric(char* fileNameOut, SemiEmpiricalModel *seModel)
936 {
937 FILE* file = NULL;
938 char buffer[1024];
939 gboolean Ok = FALSE;
940 int i;
941 int j;
942
943 file = fopen(fileNameOut, "r");
944 if(!file) return FALSE;
945 if(!fgets(buffer,BSIZE,file)) { fclose(file); return FALSE;}/* first line for energy in Hartree*/
946 if(!fgets(buffer,BSIZE,file)) { fclose(file); return FALSE;}/* dipole in au */
947 for(i=0;i<seModel->molecule.nAtoms;i++)
948 {
949 if(!fgets(buffer,BSIZE,file))break;
950 for(j=0;j<strlen(buffer);j++) if(buffer[j]=='D' || buffer[j]=='d') buffer[j] ='E';
951 if(sscanf(buffer,"%lf %lf %lf",
952 &seModel->molecule.gradient[0][i],
953 &seModel->molecule.gradient[1][i],
954 &seModel->molecule.gradient[2][i]
955 )!=3)
956 {
957 fclose(file);
958 return FALSE;
959 }
960 for(j=0;j<3;j++) seModel->molecule.gradient[j][i] *= AUTOKCAL/BOHR_TO_ANG;
961 for(j=0;j<3;j++) seModel->molecule.gradient[j][i] = - seModel->molecule.gradient[j][i];
962 }
963 Ok = TRUE;
964 fclose(file);
965 return Ok;
966 }
967 /*****************************************************************************/
runOneGeneric(SemiEmpiricalModel * seModel,char * keyWords)968 char* runOneGeneric(SemiEmpiricalModel* seModel, char* keyWords)
969 {
970 FILE* file = NULL;
971 FILE* fileSH = NULL;
972 char* fileNameIn = NULL;
973 char* fileNameOut = NULL;
974 char* fileNameSH = NULL;
975 char buffer[1024];
976 MoleculeSE* mol = &seModel->molecule;
977 char* NameCommandGeneric = seModel->method;
978 int rank = 0;
979 int type = 0;
980 #ifdef ENABLE_MPI
981 MPI_Comm_rank( MPI_COMM_WORLD,&rank);
982 #endif
983 #ifdef OS_WIN32
984 char c='%';
985 #endif
986
987 if(mol->nAtoms<1) return fileNameOut;
988 #ifndef OS_WIN32
989 fileNameSH = g_strdup_printf("%s%sGenericOne%d.sh",seModel->workDir,G_DIR_SEPARATOR_S,rank);
990 #else
991 fileNameSH = g_strdup_printf("%s%sGenericOne%d.bat",seModel->workDir,G_DIR_SEPARATOR_S,rank);
992 #endif
993 fileSH = fopen(fileNameSH, "w");
994 if(!fileSH) return FALSE;
995 #ifdef OS_WIN32
996 fprintf(fileSH,"@echo off\n");
997 #endif
998
999 fileNameIn = g_strdup_printf("%s%sOne%d.inp",seModel->workDir,G_DIR_SEPARATOR_S,rank);
1000 fileNameOut = g_strdup_printf("%s%sOne%d.out",seModel->workDir,G_DIR_SEPARATOR_S,rank);
1001
1002 file = fopen(fileNameIn, "w");
1003 if(!file)
1004 {
1005 if(fileNameIn) free(fileNameIn);
1006 if(fileNameOut) free(fileNameOut);
1007 if(fileNameSH) free(fileNameSH);
1008 return FALSE;
1009 }
1010 /*
1011 fprintf(file,"# ======================================================\n");
1012 fprintf(file,"# Generic input file made in CChemI\n");
1013 fprintf(file,"# ======================================================\n");
1014 */
1015 if(strstr(keyWords,"ENGRAD")) type = 1;
1016 fprintf(file,"%d\n",type);
1017 addMoleculeSEToFile(mol,file);
1018 fclose(file);
1019
1020 #ifndef OS_WIN32
1021 fprintf(fileSH,"%s %s %s",NameCommandGeneric,fileNameIn,fileNameOut);
1022 fclose(fileSH);
1023 sprintf(buffer,"chmod u+x %s",fileNameSH);
1024 system(buffer);
1025 system(fileNameSH);
1026 #else
1027 fprintf(fileSH,"\"%s\" \"%s\" \"%s\"",NameCommandGeneric,fileNameIn,fileNameOut);
1028 fclose(fileSH);
1029 sprintf(buffer,"\"%s\"",fileNameSH);
1030 system(buffer);
1031 #endif
1032 unlink(fileNameIn);
1033 unlink(fileNameSH);
1034 if(fileNameIn) free(fileNameIn);
1035 if(fileNameSH) free(fileNameSH);
1036 return fileNameOut;
1037 }
1038 /**********************************************************************/
newGenericModel(char * method,char * dirName,SemiEmpiricalModelConstraints constraints)1039 static SemiEmpiricalModel newGenericModel(char* method, char* dirName, SemiEmpiricalModelConstraints constraints)
1040 {
1041 /* method = nameCommand */
1042 SemiEmpiricalModel seModel = newSemiEmpiricalModel(method, dirName, constraints);
1043
1044 seModel.klass->calculateGradient = calculateGradientGeneric;
1045 seModel.klass->calculateEnergy = calculateEnergyGeneric;
1046
1047 return seModel;
1048 }
1049 /**********************************************************************/
calculateGradientGeneric(SemiEmpiricalModel * seModel)1050 static void calculateGradientGeneric(SemiEmpiricalModel* seModel)
1051 {
1052 int i;
1053 int j;
1054 MoleculeSE m = seModel->molecule;
1055 char* keyWords = NULL;
1056 char* fileOut = NULL;
1057 if(!seModel) return;
1058 if(seModel->molecule.nAtoms<1) return;
1059 if(!seModel->method) return;
1060 keyWords = g_strdup_printf("%s ENGRAD ",seModel->method);
1061 fileOut = runOneGeneric(seModel, keyWords);
1062
1063 if(fileOut)
1064 {
1065 for(j=0;j<3;j++)
1066 for( i=0; i<m.nAtoms;i++)
1067 m.gradient[j][i] = 0.0;
1068 if(!getGradientGeneric(fileOut, seModel))
1069 {
1070 #ifdef OS_WIN32
1071 char* comm = g_strdup_printf("type %s",fileOut);
1072 #else
1073 char* comm = g_strdup_printf("cat %s",fileOut);
1074 #endif
1075 printf(("Problem : I cannot caculate the Gradient... "));
1076 printf(("Calculation Stopped "));
1077 system(comm);
1078 free(fileOut);
1079 free(comm);
1080 exit(1);
1081 return;
1082 }
1083 getEnergyGeneric(fileOut, &seModel->molecule.energy);
1084 getDipoleGeneric(fileOut, seModel->molecule.dipole);
1085 free(fileOut);
1086 }
1087
1088 }
1089 /**********************************************************************/
calculateEnergyGeneric(SemiEmpiricalModel * seModel)1090 static void calculateEnergyGeneric(SemiEmpiricalModel* seModel)
1091 {
1092 char* keyWords = NULL;
1093 char* fileOut = NULL;
1094 if(!seModel) return;
1095 if(seModel->molecule.nAtoms<1) return;
1096 if(!seModel->method) return;
1097 keyWords = g_strdup_printf("%s ",seModel->method);
1098 fileOut = runOneGeneric(seModel, keyWords);
1099 if(fileOut)
1100 {
1101 getEnergyGeneric(fileOut, &seModel->molecule.energy);
1102 getDipoleGeneric(fileOut, seModel->molecule.dipole);
1103 free(fileOut);
1104 }
1105
1106 }
1107 /**********************************************************************/
createGenericModel(GeomDef * geom,gint Natoms,gint charge,gint spin,gchar * method,gchar * dirName,SemiEmpiricalModelConstraints constraints)1108 SemiEmpiricalModel createGenericModel (GeomDef* geom,gint Natoms,gint charge, gint spin, gchar* method, gchar* dirName, SemiEmpiricalModelConstraints constraints)
1109 {
1110 SemiEmpiricalModel seModel = newGenericModel(method,dirName, constraints);
1111 seModel.molecule = createMoleculeSE(geom,Natoms, charge, spin,TRUE);
1112 setRattleConstraintsParameters(&seModel);
1113
1114 return seModel;
1115 }
1116