Make parameter assignment args const char* and reduce strcpy
[namd.git] / src / Molecule.C
1 /**
2 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
3 ***  The Board of Trustees of the University of Illinois.
4 ***  All rights reserved.
5 **/
6
7 // Molecule.C is compiled twice!
8 // MOLECULE2_C undefined only compiles first half of file
9 // MOLECULE2_C defined only compiles second half of file
10 // This is shameful but it works.  Molecule needs refactoring badly.
11
12 /*
13    The class Molecule is used to hold all of the structural information
14    for a simulation.  This information is read in from a .psf file and
15    cross checked with the Parameters object passed in.  All of the structural
16    information is then stored in arrays for use.
17 */
18
19 #include "largefiles.h"  // must be first!
20
21 #include <stdio.h>
22 #include <string.h>
23 #include <stdlib.h>
24 #include <ctype.h>
25
26 #include "InfoStream.h"
27 #include "Molecule.h"
28 #include "strlib.h"
29 #include "MStream.h"
30 #include "Communicate.h"
31 #include "Node.h"
32 #include "ObjectArena.h"
33 #include "Parameters.h"
34 #include "PDB.h"
35 #include "SimParameters.h"
36 #include "Hydrogen.h"
37 #include "UniqueSetIter.h"
38 #include "charm++.h"
39 /* BEGIN gf */
40 #include "ComputeGridForce.h"
41 #include "GridForceGrid.h"
42
43 #include "MGridforceParams.h"
44 /* END gf */
45
46 #define MIN_DEBUG_LEVEL 3
47 //#define DEBUGM
48 #include "Debug.h"
49
50 #include "CompressPsf.h"
51 #include "ParallelIOMgr.h"
52 #include <deque>
53 #include <algorithm>
54
55 #ifndef M_PI
56 #define M_PI 3.14159265358979323846
57 #endif
58
59 #ifndef GROMACS_PAIR
60 #define GROMACS_PAIR 1
61 #endif
62
63 #ifndef GROMACS_EXCLUSIONS
64 #define GROMACS_EXCLUSIONS 1
65 #endif
66
67 using namespace std;
68
69 #ifndef MOLECULE2_C  // first object file
70
71 #ifdef MEM_OPT_VERSION
72 template int lookupCstPool<AtomSignature>(const vector<AtomSignature>&, const AtomSignature&);
73 template int lookupCstPool<ExclusionSignature>(const vector<ExclusionSignature>&, const ExclusionSignature&);
74 #endif
75
76 int ResidueLookupElem::lookup(
77         const char *segid, int resid, int *begin, int *end) const {
78     const ResidueLookupElem *elem = this;
79     int rval = -1;  // error
80     while ( elem && strcasecmp(elem->mySegid,segid) ) elem = elem->next;
81     if ( elem && (resid >= elem->firstResid) && (resid <= elem->lastResid) ) {
82       *begin = elem->atomIndex[resid - elem->firstResid];
83       *end = elem->atomIndex[resid - elem->firstResid + 1];
84       rval = 0;  // no error
85     }
86     return rval;
87 }
88
89 ResidueLookupElem* ResidueLookupElem::append(
90         const char *segid, int resid, int aid) {
91     ResidueLookupElem *rval = this;
92     if ( firstResid == -1 ) {  // nothing added yet
93       strcpy(mySegid,segid);
94       firstResid = resid;
95       lastResid = resid;
96       atomIndex.add(aid);
97       atomIndex.add(aid+1);
98     } else if ( ! strcasecmp(mySegid,segid) ) {  // same segid
99       if ( resid == lastResid ) {  // same resid
100         atomIndex[lastResid - firstResid + 1] = aid + 1;
101       } else if ( resid < lastResid ) {  // error
102         // We can work around this by creating a new segment.
103         iout << iWARN << "Residue " << resid <<
104           " out of order in segment " << segid <<
105           ", lookup for additional residues in this segment disabled.\n" << endi;
106         rval = next = new ResidueLookupElem;
107         next->append(segid,resid,aid);
108       } else {  // new resid
109         for ( ; lastResid < resid; ++lastResid ) atomIndex.add(aid);
110         atomIndex[lastResid - firstResid + 1] = aid + 1;
111       }
112     } else {  // new segid
113       rval = next = new ResidueLookupElem;
114       next->append(segid,resid,aid);
115     }
116     return rval;
117 }
118
119
120 //  Lookup atom id from segment, residue, and name
121 int Molecule::get_atom_from_name(
122         const char *segid, int resid, const char *aname) const {
123
124   if (atomNames == NULL || resLookup == NULL)
125   {
126     NAMD_die("Tried to find atom from name on node other than node 0");
127   }
128
129   int i = 0;
130   int end = 0;
131   if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
132   for ( ; i < end; ++i ) {
133     #ifdef MEM_OPT_VERSION    
134     Index idx = atomNames[i].atomnameIdx;
135     if(!strcasecmp(aname, atomNamePool[idx])) return i;
136     #else
137     if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i;
138     #endif
139   }
140   return -1;
141 }
142
143 //  Lookup number of atoms in residue from segment and residue
144 int Molecule::get_residue_size(
145         const char *segid, int resid) const {
146
147   if (atomNames == NULL || resLookup == NULL)
148   {
149     NAMD_die("Tried to find atom from name on node other than node 0");
150   }
151   int i = 0;
152   int end = 0;
153   if ( resLookup->lookup(segid,resid,&i,&end) ) return 0;
154   return ( end - i );
155 }
156
157 //  Lookup atom id from segment, residue, and index in residue
158 int Molecule::get_atom_from_index_in_residue(
159         const char *segid, int resid, int index) const {
160
161   if (atomNames == NULL || resLookup == NULL)
162   {
163     NAMD_die("Tried to find atom from name on node other than node 0");
164   }
165   int i = 0;
166   int end = 0;
167   if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
168   if ( index >= 0 && index < ( end - i ) ) return ( index + i );
169   return -1;
170 }
171
172 /************************************************************************/
173 /*                  */
174 /*      FUNCTION initialize  */
175 /*                  */
176 /*  This is the initializer for the Molecule class.  It simply sets */
177 /*  the counts for all the various parameters to 0 and sets the pointers*/
178 /*  to the arrays that will store these parameters to NULL, since they  */
179 /*  have not been allocated yet.          */
180 /*                  */
181 /************************************************************************/
182
183 void Molecule::initialize(SimParameters *simParams, Parameters *param)
184 {
185   if ( sizeof(int32) != 4 ) { NAMD_bug("sizeof(int32) != 4"); }
186   this->simParams = simParams;
187   this->params = param;
188
189   /*  Initialize array pointers to NULL  */
190   atoms=NULL;
191   atomNames=NULL;
192   resLookup=NULL;
193
194   // DRUDE
195   is_lonepairs_psf = 1;  // anticipate having lone pair hosts
196   is_drude_psf = 0;  // assume not Drude model
197   drudeConsts=NULL;
198   lphosts=NULL;  // might have lone pair hosts without Drude!
199   anisos=NULL;
200   tholes=NULL;
201   lphostIndexes=NULL;
202   // DRUDE
203
204   //LCPO
205   lcpoParamType = NULL;
206
207   //for compressing molecule info
208   atomSegResids=NULL;
209
210   if ( simParams->globalForcesOn ) {
211     resLookup = new ResidueLookupElem;
212   }
213
214   #ifdef MEM_OPT_VERSION
215   eachAtomSig = NULL;
216   atomSigPoolSize = 0;
217   atomSigPool = NULL;
218   massPoolSize = 0;
219   atomMassPool = NULL;
220   eachAtomMass = NULL;
221   chargePoolSize = 0;
222   atomChargePool = NULL;
223   eachAtomCharge = NULL;
224   #else
225   bonds=NULL;
226   angles=NULL;
227   dihedrals=NULL;
228   impropers=NULL;
229   crossterms=NULL;
230   #endif
231
232   donors=NULL;
233   acceptors=NULL;
234   
235
236   #ifndef MEM_OPT_VERSION      
237   tmpArena=NULL;
238   exclusions=NULL;
239   bondsWithAtom=NULL;
240   bondsByAtom=NULL;
241   anglesByAtom=NULL;
242   dihedralsByAtom=NULL;
243   impropersByAtom=NULL;
244   crosstermsByAtom=NULL;
245   // JLai
246   gromacsPairByAtom=NULL;
247   // End of JLai
248   // DRUDE
249   tholesByAtom=NULL;
250   anisosByAtom=NULL;
251   // DRUDE
252   #endif
253
254   #ifdef MEM_OPT_VERSION
255   exclSigPool = NULL;
256   exclChkSigPool = NULL;
257   exclSigPoolSize = 0;
258   eachAtomExclSig = NULL;
259
260   fixedAtomsSet = NULL;
261   constrainedAtomsSet = NULL;
262   #else
263   exclusionsByAtom=NULL;
264   fullExclusionsByAtom=NULL;
265   modExclusionsByAtom=NULL;
266   all_exclusions=NULL;
267   #endif
268
269   langevinParams=NULL;
270   fixedAtomFlags=NULL;
271
272   #ifdef MEM_OPT_VERSION  
273   clusterSigs=NULL;
274   #else
275   cluster=NULL;  
276   #endif
277   clusterSize=NULL;
278
279   exPressureAtomFlags=NULL;
280   rigidBondLengths=NULL;
281   consIndexes=NULL;
282   consParams=NULL;
283   /* BEGIN gf */
284   gridfrcIndexes=NULL;
285   gridfrcParams=NULL;
286   gridfrcGrid=NULL;
287   numGridforces=NULL;
288   /* END gf */
289   stirIndexes=NULL;
290   stirParams=NULL;
291   movDragIndexes=NULL;
292   movDragParams=NULL;
293   rotDragIndexes=NULL;
294   rotDragParams=NULL;
295   consTorqueIndexes=NULL;
296   consTorqueParams=NULL;
297   consForceIndexes=NULL;
298   consForce=NULL;
299 //fepb
300   fepAtomFlags=NULL;
301 //fepe
302 //soluteScaling
303   ssAtomFlags=NULL;
304   ss_vdw_type=NULL;
305   ss_index=NULL;
306 //soluteScaling
307   nameArena = new ObjectArena<char>;
308   // nameArena->setAlignment(8);
309   // arena->setAlignment(32);
310   #ifndef MEM_OPT_VERSION
311   arena = new ObjectArena<int32>;
312   exclArena = new ObjectArena<char>;
313   #endif
314   // exclArena->setAlignment(32);
315
316   /*  Initialize counts to 0 */
317   numAtoms=0;
318   numRealBonds=0;
319   numBonds=0;
320   numAngles=0;
321   numDihedrals=0;
322   numImpropers=0;
323   numTholes=0;
324   numAnisos=0;
325   numCrossterms=0;
326   // JLai
327   numLJPair=0;
328   // End of JLai
329   numDonors=0;
330   numAcceptors=0;
331   numExclusions=0;
332
333   // DRUDE
334   numLonepairs=0;
335   numDrudeAtoms=0;
336   numLphosts=0;
337   numAnisos=0;
338   // DRUDE
339
340   numConstraints=0;
341   numStirredAtoms=0;
342   numMovDrag=0;
343   numRotDrag=0;
344   numConsTorque=0;
345   numConsForce=0;
346   numFixedAtoms=0;
347   numFixedGroups=0;
348   numExPressureAtoms=0;
349   numRigidBonds=0;
350   numFixedRigidBonds=0;
351   numMultipleDihedrals=0;
352   numMultipleImpropers=0;
353   numCalcBonds=0;
354   numCalcAngles=0;
355   numCalcDihedrals=0;
356   numCalcImpropers=0;
357   numCalcTholes=0;
358   numCalcAnisos=0;
359   numCalcCrossterms=0;
360   numCalcExclusions=0;
361   numCalcFullExclusions=0;
362   // JLai
363   numCalcLJPair=0;
364   // End of JLai
365
366 //fepb
367   numFepInitial = 0;
368   numFepFinal = 0;
369 //fepe
370
371   //fields related with pluginIO-based loading molecule structure
372   occupancy = NULL;
373   bfactor = NULL;
374
375   qmElementArray=0;
376   qmDummyElement=0;
377   qmGrpSizes=0;
378   qmAtomGroup=0;
379   qmAtmChrg=0;
380   qmAtmIndx=0;
381   qmGrpID=0;
382   qmGrpChrg=0;
383   qmGrpMult=0;
384   qmGrpNumBonds=0;
385   qmMMBond=0;
386   qmGrpBonds=0;
387   qmMMBondedIndx=0;
388   qmMMChargeTarget=0;
389   qmMMNumTargs=0;
390   qmDummyBondVal=0;
391   qmMeMMindx=0;
392   qmMeQMGrp=0;
393   qmCustomPCIdxs=0;
394   qmCustPCSizes=0;
395   qmLSSSize=0;
396   qmLSSIdxs=0;
397   qmLSSMass=0;
398   qmLSSRefIDs=0;
399   qmLSSRefMass=0;
400   qmLSSRefSize=0;
401   qmNumBonds=0;
402   cSMDindex=0;
403   cSMDindxLen=0;
404   cSMDpairs=0;
405   cSMDKs=0;
406   cSMDVels=0;
407   cSMDcoffs=0;
408   cSMDnumInst=0;
409   
410   goInit();
411 }
412
413 /*      END OF FUNCTION initialize */
414
415 /************************************************************************/
416 /*                  */
417 /*      FUNCTION Molecule        */
418 /*                  */
419 /*  This is the constructor for the Molecule class. */
420 /*                  */
421 /************************************************************************/
422
423 Molecule::Molecule(SimParameters *simParams, Parameters *param)
424 {
425   initialize(simParams,param);
426 }
427
428 /************************************************************************/
429 /*                  */
430 /*      FUNCTION Molecule        */
431 /*                  */
432 /*  This is the constructor for the Molecule class from CHARMM/XPLOR files. */
433 /*                  */
434 /************************************************************************/
435
436 Molecule::Molecule(SimParameters *simParams, Parameters *param, char *filename, ConfigList *cfgList)
437 {
438   initialize(simParams,param);
439
440 #ifdef MEM_OPT_VERSION
441   if(simParams->useCompressedPsf)
442       read_mol_signatures(filename, param, cfgList);
443 #else
444         read_psf_file(filename, param); 
445  //LCPO
446   if (simParams->LCPOOn)
447     assignLCPOTypes( 0 );
448 #endif      
449  }
450
451 /************************************************************************/
452 /*                                                                      */
453 /*      FUNCTION Molecule                                               */
454 /*                                                                      */
455 /*  This is the constructor for the Molecule class from plugin IO.      */
456 /*                                                                      */
457 /************************************************************************/
458 Molecule::Molecule(SimParameters *simParams, Parameters *param, molfile_plugin_t *pIOHdl, void *pIOFileHdl, int natoms)
459 {
460 #ifdef MEM_OPT_VERSION
461   NAMD_die("Sorry, plugin IO is not supported in the memory optimized version.");
462 #else
463     initialize(simParams, param);
464     numAtoms = natoms;
465     int optflags = MOLFILE_BADOPTIONS;
466     molfile_atom_t *atomarray = (molfile_atom_t *) malloc(natoms*sizeof(molfile_atom_t));
467     memset(atomarray, 0, natoms*sizeof(molfile_atom_t));
468
469     //1a. read basic atoms information
470     int rc = pIOHdl->read_structure(pIOFileHdl, &optflags, atomarray);
471     if (rc != MOLFILE_SUCCESS && rc != MOLFILE_NOSTRUCTUREDATA) {
472         free(atomarray);
473         NAMD_die("ERROR: plugin failed reading structure data");
474     }
475     if(optflags == MOLFILE_BADOPTIONS) {
476         free(atomarray);
477         NAMD_die("ERROR: plugin didn't initialize optional data flags");
478     }
479     if(optflags & MOLFILE_OCCUPANCY) {
480         setOccupancyData(atomarray);
481     }
482     if(optflags & MOLFILE_BFACTOR) {
483         setBFactorData(atomarray);
484     }
485     //1b. load basic atoms information to the molecule object
486     plgLoadAtomBasics(atomarray);    
487     free(atomarray);
488
489     //2a. read bonds
490     //indices are one-based in read_bonds
491     int *from, *to;
492     float *bondorder;
493     int *bondtype, nbondtypes;
494     char **bondtypename;
495     if(pIOHdl->read_bonds!=NULL) {
496         if(pIOHdl->read_bonds(pIOFileHdl, &numBonds, &from, &to, &bondorder,
497                                  &bondtype, &nbondtypes, &bondtypename)){
498             NAMD_die("ERROR: failed reading bond information.");
499         }
500     }    
501     //2b. load bonds information to the molecule object
502     if(numBonds!=0) {
503         plgLoadBonds(from,to);
504     }
505
506     //3a. read other bonded structures
507     int *plgAngles, *plgDihedrals, *plgImpropers, *plgCterms;
508     int ctermcols, ctermrows;
509     int *angletypes, numangletypes, *dihedraltypes, numdihedraltypes;
510     int *impropertypes, numimpropertypes; 
511     char **angletypenames, **dihedraltypenames, **impropertypenames;
512
513     plgAngles=plgDihedrals=plgImpropers=plgCterms=NULL;
514     if(pIOHdl->read_angles!=NULL) {
515         if(pIOHdl->read_angles(pIOFileHdl,
516                   &numAngles, &plgAngles,
517                   &angletypes, &numangletypes, &angletypenames,
518                   &numDihedrals, &plgDihedrals,
519                   &dihedraltypes, &numdihedraltypes, &dihedraltypenames,
520                   &numImpropers, &plgImpropers,
521                   &impropertypes, &numimpropertypes, &impropertypenames,
522                   &numCrossterms, &plgCterms, &ctermcols, &ctermrows)) {
523             NAMD_die("ERROR: failed reading angle information.");
524         }
525     }
526     //3b. load other bonded structures to the molecule object
527     if(numAngles!=0) plgLoadAngles(plgAngles);
528     if(numDihedrals!=0) plgLoadDihedrals(plgDihedrals);
529     if(numImpropers!=0) plgLoadImpropers(plgImpropers);
530     if(numCrossterms!=0) plgLoadCrossterms(plgCterms);
531
532   numRealBonds = numBonds;
533   build_atom_status();
534   //LCPO
535   if (simParams->LCPOOn)
536     assignLCPOTypes( 2 );
537 #endif
538 }
539
540 /*      END OF FUNCTION Molecule      */
541
542 /************************************************************************/
543 /*                  */
544 /*        FUNCTION Molecule      */
545 /*                  */
546 /*  This is the destructor for the class Molecule.  It simply frees */
547 /*  the memory allocated for each of the arrays used to store the       */
548 /*  structure information.            */
549 /*                  */
550 /************************************************************************/
551
552 Molecule::~Molecule()
553 {
554   /*  Check to see if each array was ever allocated.  If it was   */
555   /*  then free it            */
556   if (atoms != NULL)
557     delete [] atoms;
558
559   if (atomNames != NULL)
560   {
561     // subarrarys allocated from arena - automatically deleted
562     delete [] atomNames;
563   }
564   delete nameArena;
565
566   if (resLookup != NULL)
567     delete resLookup;
568
569   // DRUDE: free arrays read from PSF
570   if (drudeConsts != NULL) delete [] drudeConsts;
571   if (lphosts != NULL) delete [] lphosts;
572   if (anisos != NULL) delete [] anisos;
573   if (tholes != NULL) delete [] tholes;
574   if (lphostIndexes != NULL) delete [] lphostIndexes;
575   // DRUDE
576
577   //LCPO
578   if (lcpoParamType != NULL) delete [] lcpoParamType;
579
580   #ifdef MEM_OPT_VERSION
581   if(eachAtomSig) delete [] eachAtomSig;
582   if(atomSigPool) delete [] atomSigPool;
583   #else
584   if (bonds != NULL)
585     delete [] bonds;
586
587   if (angles != NULL)
588     delete [] angles;
589
590   if (dihedrals != NULL)
591     delete [] dihedrals;
592
593   if (impropers != NULL)
594     delete [] impropers;
595
596   if (crossterms != NULL)
597     delete [] crossterms;
598
599   if (exclusions != NULL)
600     delete [] exclusions;
601   #endif
602
603   if (donors != NULL)
604     delete [] donors;
605
606   if (acceptors != NULL)
607     delete [] acceptors;  
608
609   #ifdef MEM_OPT_VERSION
610   if(exclSigPool) delete [] exclSigPool;
611   if(exclChkSigPool) delete [] exclChkSigPool;
612   if(eachAtomExclSig) delete [] eachAtomExclSig;
613   if(fixedAtomsSet) delete fixedAtomsSet;
614   if(constrainedAtomsSet) delete constrainedAtomsSet;
615   #else
616   if (bondsByAtom != NULL)
617        delete [] bondsByAtom;
618   
619   if (anglesByAtom != NULL)
620        delete [] anglesByAtom;
621   
622   if (dihedralsByAtom != NULL)
623        delete [] dihedralsByAtom;
624   
625   if (impropersByAtom != NULL)
626        delete [] impropersByAtom;
627   
628   if (crosstermsByAtom != NULL)
629        delete [] crosstermsByAtom;  
630
631   if (exclusionsByAtom != NULL)
632        delete [] exclusionsByAtom;
633   
634   if (fullExclusionsByAtom != NULL)
635        delete [] fullExclusionsByAtom;
636   
637   if (modExclusionsByAtom != NULL)
638        delete [] modExclusionsByAtom;
639   
640   if (all_exclusions != NULL)
641        delete [] all_exclusions;
642
643   // JLai
644   if (gromacsPairByAtom != NULL)
645       delete [] gromacsPairByAtom;
646   // End of JLai
647
648   // DRUDE
649   if (tholesByAtom != NULL)
650        delete [] tholesByAtom;
651   if (anisosByAtom != NULL)
652        delete [] anisosByAtom;
653   // DRUDE
654   #endif
655
656   //LCPO
657   if (lcpoParamType != NULL)
658     delete [] lcpoParamType;
659
660   if (fixedAtomFlags != NULL)
661        delete [] fixedAtomFlags;
662
663   if (stirIndexes != NULL)
664     delete [] stirIndexes;
665
666
667   #ifdef MEM_OPT_VERSION
668   if(clusterSigs != NULL){      
669       delete [] clusterSigs;
670   }  
671   #else
672   if (cluster != NULL)
673        delete [] cluster;  
674   #endif
675   if (clusterSize != NULL)
676        delete [] clusterSize;
677
678   if (exPressureAtomFlags != NULL)
679        delete [] exPressureAtomFlags;
680
681   if (rigidBondLengths != NULL)
682        delete [] rigidBondLengths;
683
684 //fepb
685   if (fepAtomFlags != NULL)
686        delete [] fepAtomFlags;
687 //fepe
688
689 //soluteScaling
690   if (ssAtomFlags != NULL)
691        delete [] ssAtomFlags;
692   if (ss_vdw_type != NULL)
693        delete [] ss_vdw_type;
694   if (ss_index != NULL)
695        delete [] ss_index;
696 //soluteScaling
697
698   if (qmAtomGroup != NULL)
699        delete [] qmAtomGroup;
700   
701   if (qmAtmIndx != NULL)
702        delete [] qmAtmIndx;
703   
704   if (qmAtmChrg != NULL)
705        delete [] qmAtmChrg;
706   
707   
708   if (qmGrpNumBonds != NULL)
709        delete [] qmGrpNumBonds;
710   
711   if (qmGrpSizes != NULL)
712        delete [] qmGrpSizes;
713   
714   if (qmDummyBondVal != NULL)
715        delete [] qmDummyBondVal;
716   
717   if (qmMMNumTargs != NULL)
718        delete [] qmMMNumTargs;
719   
720   if (qmGrpID != NULL)
721        delete [] qmGrpID;
722   
723   if (qmGrpChrg != NULL)
724        delete [] qmGrpChrg;
725   
726   if (qmGrpMult != NULL)
727        delete [] qmGrpMult;
728   
729   if (qmMeMMindx != NULL)
730        delete [] qmMeMMindx;
731   
732   if (qmMeQMGrp != NULL)
733        delete [] qmMeQMGrp;
734   
735   if (qmLSSSize != NULL)
736        delete [] qmLSSSize;
737   
738   if (qmLSSIdxs != NULL)
739        delete [] qmLSSIdxs;
740   
741   if (qmLSSMass != NULL)
742        delete [] qmLSSMass;
743   
744   if (qmLSSRefSize != NULL)
745        delete [] qmLSSRefSize;
746   
747   if (qmLSSRefIDs != NULL)
748        delete [] qmLSSRefIDs;
749   
750   if (qmLSSRefMass != NULL)
751        delete [] qmLSSRefMass;
752   
753   if (qmMMBond != NULL) {
754       for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
755           if (qmMMBond[grpIndx] != NULL)
756               delete [] qmMMBond[grpIndx];
757       }
758       delete [] qmMMBond;
759   }
760   
761   if (qmGrpBonds != NULL) {
762       for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
763           if (qmGrpBonds[grpIndx] != NULL)
764               delete [] qmGrpBonds[grpIndx];
765       }
766       delete [] qmGrpBonds;
767   }
768   
769   if (qmMMBondedIndx != NULL) {
770       for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
771           if (qmMMBondedIndx[grpIndx] != NULL)
772               delete [] qmMMBondedIndx[grpIndx];
773       }
774       delete [] qmMMBondedIndx;
775   }
776   
777   if (qmMMChargeTarget != NULL) {
778       for(int grpIndx = 0 ; grpIndx < qmNumBonds; grpIndx++) {
779           if (qmMMChargeTarget[grpIndx] != NULL)
780               delete [] qmMMChargeTarget[grpIndx];
781       }
782       delete [] qmMMChargeTarget;
783   }
784   
785     if (qmElementArray != NULL){
786         for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
787             if (qmElementArray[atmInd] != NULL)
788                 delete [] qmElementArray[atmInd];
789         }
790         delete [] qmElementArray;
791     }
792     
793     if (qmDummyElement != NULL){
794         for(int atmInd = 0 ; atmInd < numAtoms; atmInd++) {
795             if (qmDummyElement[atmInd] != NULL)
796                 delete [] qmDummyElement[atmInd];
797         }
798         delete [] qmDummyElement;
799     }
800
801     if (qmCustPCSizes != NULL){
802         delete [] qmCustPCSizes;
803     }
804     
805     if (qmCustomPCIdxs != NULL){
806         delete [] qmCustomPCIdxs;
807     }
808     
809     if (cSMDindex != NULL) {
810       for(int grpIndx = 0 ; grpIndx < qmNumGrps; grpIndx++) {
811           if (cSMDindex[grpIndx] != NULL)
812               delete [] cSMDindex[grpIndx];
813       }
814       delete [] cSMDindex;
815     }
816     
817     if (cSMDpairs != NULL) {
818       for(int instIndx = 0 ; instIndx < cSMDnumInst; instIndx++) {
819           if (cSMDpairs[instIndx] != NULL)
820               delete [] cSMDpairs[instIndx];
821       }
822       delete [] cSMDpairs;
823     }
824     
825     if (cSMDindxLen != NULL)
826         delete [] cSMDindxLen;
827     if (cSMDKs != NULL)
828         delete [] cSMDKs;
829     if (cSMDVels != NULL)
830         delete [] cSMDVels;
831     if (cSMDcoffs != NULL)
832         delete [] cSMDcoffs;
833     
834   #ifndef MEM_OPT_VERSION
835   delete arena;
836   delete exclArena;
837   #endif
838 }
839 /*      END OF FUNCTION Molecule      */
840
841 #ifndef MEM_OPT_VERSION
842
843 //===Non-memory optimized version of functions that read Molecule file===//
844
845 /************************************************************************/
846 /*                  */
847 /*        FUNCTION read_psf_file      */
848 /*                  */
849 /*   INPUTS:                */
850 /*  fname - Name of the .psf file to read        */
851 /*  params - pointer to Parameters object to use to obtain          */
852 /*     parameters for vdWs, bonds, etc.      */
853 /*                  */
854 /*  This function reads a .psf file in.  This is where just about   */
855 /*   all of the structural information for this class comes from.  The  */
856 /*   .psf file contains descriptions of the atom, bonds, angles,        */
857 /*   dihedrals, impropers, and exclusions.  The parameter object is     */
858 /*   used to look up parameters for each of these entities.    */
859 /*                  */
860 /************************************************************************/
861
862 void Molecule::read_psf_file(char *fname, Parameters *params)
863 {
864   char err_msg[512];  //  Error message for NAMD_die
865   char buffer[512];  //  Buffer for file reading
866   int i;      //  Loop counter
867   int NumTitle;    //  Number of Title lines in .psf file
868   FILE *psf_file;    //  pointer to .psf file
869   int ret_code;    //  ret_code from NAMD_read_line calls
870
871   /* Try and open the .psf file           */
872   if ( (psf_file = Fopen(fname, "r")) == NULL)
873   {
874     sprintf(err_msg, "UNABLE TO OPEN .psf FILE %s", fname);
875     NAMD_die(err_msg);
876   }
877
878   /*  Read till we have the first non-blank line of file    */
879   ret_code = NAMD_read_line(psf_file, buffer);
880
881   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
882   {
883     ret_code = NAMD_read_line(psf_file, buffer);
884   }
885
886   /*  Check to see if we dropped out of the loop because of a     */
887   /*  read error.  This shouldn't happen unless the file is empty */
888   if (ret_code!=0)
889   {
890     sprintf(err_msg, "EMPTY .psf FILE %s", fname);
891     NAMD_die(err_msg);
892   }
893
894   /*  The first non-blank line should contain the word "psf".    */
895   /*  If we can't find it, die.               */
896   if (!NAMD_find_word(buffer, "psf"))
897   {
898     sprintf(err_msg, "UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
899        fname);
900     NAMD_die(err_msg);
901   }
902
903   // DRUDE: set flag if we discover Drude PSF
904   if (NAMD_find_word(buffer, "drude"))
905   {
906     if ( ! simParams->drudeOn ) {
907       iout << iWARN << "Reading PSF supporting DRUDE without "
908         "enabling the Drude model in the simulation config file\n" << endi;
909     }
910     is_drude_psf = 1;
911   }
912   // DRUDE
913
914   /*  Read until we find the next non-blank line      */
915   ret_code = NAMD_read_line(psf_file, buffer);
916
917   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
918   {
919     ret_code = NAMD_read_line(psf_file, buffer);
920   }
921
922   /*  Check to see if we dropped out of the loop because of a     */
923   /*  read error.  This shouldn't happen unless there is nothing  */
924   /*  but the PSF line in the file        */
925   if (ret_code!=0)
926   {
927     sprintf(err_msg, "MISSING EVERYTHING BUT PSF FROM %s", fname);
928     NAMD_die(err_msg);
929   }
930
931   /*  This line should have the word "NTITLE" in it specifying    */
932   /*  how many title lines there are        */
933   if (!NAMD_find_word(buffer, "NTITLE"))
934   {
935     sprintf(err_msg,"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
936        fname);
937     NAMD_die(err_msg);
938   }
939
940   sscanf(buffer, "%d", &NumTitle);
941
942   /*  Now skip the next NTITLE non-blank lines and then read in the*/
943   /*  line which should contain NATOM        */
944   i=0;
945
946   while ( ((ret_code=NAMD_read_line(psf_file, buffer)) == 0) && 
947     (i<NumTitle) )
948   {
949     if (!NAMD_blank_string(buffer))
950       i++;
951   }
952
953   /*  Make sure we didn't exit because of a read error    */
954   if (ret_code!=0)
955   {
956     sprintf(err_msg, "FOUND EOF INSTEAD OF NATOM IN PSF FILE %s", 
957        fname);
958     NAMD_die(err_msg);
959   }
960
961   while (NAMD_blank_string(buffer))
962   {
963     NAMD_read_line(psf_file, buffer);
964   }
965
966   /*  Check to make sure we have the line we want      */
967   if (!NAMD_find_word(buffer, "NATOM"))
968   {
969     sprintf(err_msg, "DIDN'T FIND \"NATOM\" IN PSF FILE %s",
970        fname);
971     NAMD_die(err_msg);
972   }
973
974   /*  Read in the number of atoms, and then the atoms themselves  */
975   sscanf(buffer, "%d", &numAtoms);
976
977   read_atoms(psf_file, params);
978
979   /*  Read until we find the next non-blank line      */
980   ret_code = NAMD_read_line(psf_file, buffer);
981
982   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
983   {
984     ret_code = NAMD_read_line(psf_file, buffer);
985   }
986
987   /*  Check to make sure we didn't hit the EOF      */
988   if (ret_code != 0)
989   {
990     NAMD_die("EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
991   }
992
993   /*  Look for the string "NBOND"          */
994   if (!NAMD_find_word(buffer, "NBOND"))
995   {
996     NAMD_die("DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
997   }
998
999   /*  Read in the number of bonds and then the bonds themselves  */
1000   sscanf(buffer, "%d", &numBonds);
1001
1002   if (numBonds)
1003     read_bonds(psf_file, params);
1004
1005   /*  Read until we find the next non-blank line      */
1006   ret_code = NAMD_read_line(psf_file, buffer);
1007
1008   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
1009   {
1010     ret_code = NAMD_read_line(psf_file, buffer);
1011   }
1012
1013   /*  Check to make sure we didn't hit the EOF      */
1014   if (ret_code != 0)
1015   {
1016     NAMD_die("EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
1017   }
1018
1019   /*  Look for the string "NTHETA"        */
1020   if (!NAMD_find_word(buffer, "NTHETA"))
1021   {
1022     NAMD_die("DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
1023   }
1024
1025   /*  Read in the number of angles and then the angles themselves */
1026   sscanf(buffer, "%d", &numAngles);
1027
1028   if (numAngles)
1029     read_angles(psf_file, params);
1030
1031   /*  Read until we find the next non-blank line      */
1032   ret_code = NAMD_read_line(psf_file, buffer);
1033
1034   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
1035   {
1036     ret_code = NAMD_read_line(psf_file, buffer);
1037   }
1038
1039   /*  Check to make sure we didn't hit the EOF      */
1040   if (ret_code != 0)
1041   {
1042     NAMD_die("EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
1043   }
1044
1045   /*  Look for the string "NPHI"          */
1046   if (!NAMD_find_word(buffer, "NPHI"))
1047   {
1048     NAMD_die("DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
1049   }
1050
1051   /*  Read in the number of dihedrals and then the dihedrals      */
1052   sscanf(buffer, "%d", &numDihedrals);
1053
1054   if (numDihedrals)
1055     read_dihedrals(psf_file, params);
1056
1057   /*  Read until we find the next non-blank line      */
1058   ret_code = NAMD_read_line(psf_file, buffer);
1059
1060   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
1061   {
1062     ret_code = NAMD_read_line(psf_file, buffer);
1063   }
1064
1065   /*  Check to make sure we didn't hit the EOF      */
1066   if (ret_code != 0)
1067   {
1068     NAMD_die("EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
1069   }
1070
1071   /*  Look for the string "NIMPHI"        */
1072   if (!NAMD_find_word(buffer, "NIMPHI"))
1073   {
1074     NAMD_die("DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
1075   }
1076
1077   /*  Read in the number of Impropers and then the impropers  */
1078   sscanf(buffer, "%d", &numImpropers);
1079
1080   if (numImpropers)
1081     read_impropers(psf_file, params);
1082
1083   /*  Read until we find the next non-blank line      */
1084   ret_code = NAMD_read_line(psf_file, buffer);
1085
1086   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
1087   {
1088     ret_code = NAMD_read_line(psf_file, buffer);
1089   }
1090
1091   /*  Check to make sure we didn't hit the EOF      */
1092   if (ret_code != 0)
1093   {
1094     NAMD_die("EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
1095   }
1096
1097   /*  Look for the string "NDON"        */
1098   if (!NAMD_find_word(buffer, "NDON"))
1099   {
1100     NAMD_die("DID NOT FIND NDON AFTER ATOM LIST IN PSF");
1101   }
1102
1103   /*  Read in the number of hydrogen bond donors and then the donors */
1104   sscanf(buffer, "%d", &numDonors);
1105
1106   if (numDonors)
1107     read_donors(psf_file);
1108
1109   /*  Read until we find the next non-blank line      */
1110   ret_code = NAMD_read_line(psf_file, buffer);
1111
1112   while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
1113   {
1114     ret_code = NAMD_read_line(psf_file, buffer);
1115   }
1116
1117   /*  Check to make sure we didn't hit the EOF      */
1118   if (ret_code != 0)
1119   {
1120     NAMD_die("EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
1121   }
1122
1123   /*  Look for the string "NACC"        */
1124   if (!NAMD_find_word(buffer, "NACC"))
1125   {
1126     NAMD_die("DID NOT FIND NACC AFTER ATOM LIST IN PSF");
1127   }
1128
1129   /*  Read in the number of hydrogen bond donors and then the donors */
1130   sscanf(buffer, "%d", &numAcceptors);
1131
1132   if (numAcceptors)
1133     read_acceptors(psf_file);
1134
1135   /*  look for the explicit non-bonded exclusion section.     */
1136   while (!NAMD_find_word(buffer, "NNB"))
1137   {
1138     ret_code = NAMD_read_line(psf_file, buffer);
1139
1140     if (ret_code != 0)
1141     {
1142       NAMD_die("EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
1143     }
1144   }
1145
1146   /*  Read in the number of exclusions and then the exclusions    */
1147   sscanf(buffer, "%d", &numExclusions);
1148
1149   if (numExclusions)
1150     read_exclusions(psf_file);
1151
1152   //
1153   // The remaining sections that we can read might be optional to the PSF.
1154   //
1155   // Possibilities are:
1156   // - Drude simulation:
1157   //   + NUMLP (probably, but necessarily?)
1158   //   + NUMANISO (required)
1159   //   + NCRTERM (optional)
1160   // - non-Drude simulation
1161   //   + NUMLP (optional)
1162   //   + NCRTERM (optional)
1163   //
1164   // Also, there might be unrecognized PSF sections that we should skip.
1165   //
1166
1167   // Keep reading until we recognize another section of PSF or find EOF.
1168   int is_found = 0;
1169   int is_found_numlp = 0;
1170   int is_found_numaniso = 0;
1171   int is_found_ncrterm = 0;
1172   while ( ! is_found ) {
1173     // Read until we find the next non-blank line
1174     do {
1175       ret_code = NAMD_read_line(psf_file, buffer);
1176     } while (ret_code == 0 && NAMD_blank_string(buffer) != 0 );
1177     // we either found EOF or we will try to match word in buffer
1178     if (ret_code != 0) {
1179       is_found = -1;  // found end of file
1180     }
1181     else if ( (is_found_numlp = NAMD_find_word(buffer, "NUMLP")) > 0) {
1182       is_found = is_found_numlp;
1183     }
1184     else if ( (is_found_numaniso = NAMD_find_word(buffer, "NUMANISO")) > 0) {
1185       is_found = is_found_numaniso;
1186     }
1187     else if ( (is_found_ncrterm = NAMD_find_word(buffer, "NCRTERM")) > 0) {
1188       is_found = is_found_ncrterm;
1189     }
1190   }
1191
1192   if (is_found_numlp) {
1193     // We found lone pair hosts.
1194     // Read the number and then read the lone pair hosts.
1195     sscanf(buffer, "%d", &numLphosts);
1196     if (numLphosts > 0) read_lphosts(psf_file);
1197
1198     // Keep reading for next keyword.
1199     is_found = 0;
1200     is_found_numaniso = 0;
1201     is_found_ncrterm = 0;
1202     while ( ! is_found ) {
1203       // Read until we find the next non-blank line
1204       do {
1205         ret_code = NAMD_read_line(psf_file, buffer);
1206       } while (ret_code == 0 && NAMD_blank_string(buffer) != 0 );
1207       // we either found EOF or we will try to match word in buffer
1208       if (ret_code != 0) {
1209         is_found = -1;  // found end of file
1210       }
1211       else if ( (is_found_numaniso = NAMD_find_word(buffer, "NUMANISO")) > 0) {
1212         is_found = is_found_numaniso;
1213       }
1214       else if ( (is_found_ncrterm = NAMD_find_word(buffer, "NCRTERM")) > 0) {
1215         is_found = is_found_ncrterm;
1216       }
1217     }
1218   }
1219
1220   if (numLphosts == 0) {
1221     // We either had no NUMLP section or we did and zero were listed.
1222     // Nevertheless, we have no lone pair hosts
1223     // so reset the simparams flag before simulation
1224     simParams->lonepairs = FALSE;
1225     is_lonepairs_psf = 0;
1226   }
1227   else if (simParams->lonepairs == FALSE /* but numLphosts > 0 */) {
1228     // Config file "lonepairs" option is (now) enabled by default.
1229     // Bad things will happen when lone pair hosts exist but "lonepairs"
1230     // in simparams is disabled. In this event, terminate with an error.
1231     NAMD_die("FOUND LONE PAIR HOSTS IN PSF WITH \"LONEPAIRS\" DISABLED IN CONFIG FILE");
1232   }
1233
1234   if (is_found_numaniso && is_drude_psf) {
1235     // We are reading a Drude PSF and found the anisotropic terms.
1236     // Read the number and then read those terms.
1237     sscanf(buffer, "%d", &numAnisos);
1238     if (numAnisos > 0) read_anisos(psf_file);
1239
1240     // Keep reading for next keyword.
1241     is_found = 0;
1242     is_found_ncrterm = 0;
1243     while ( ! is_found ) {
1244       // Read until we find the next non-blank line
1245       do {
1246         ret_code = NAMD_read_line(psf_file, buffer);
1247       } while (ret_code == 0 && NAMD_blank_string(buffer) != 0 );
1248       // we either found EOF or we will try to match word in buffer
1249       if (ret_code != 0) {
1250         is_found = -1;  // found end of file
1251       }
1252       else if ( (is_found_ncrterm = NAMD_find_word(buffer, "NCRTERM")) > 0) {
1253         is_found = is_found_ncrterm;
1254       }
1255     }
1256   }
1257   else if (is_drude_psf /* but not is_found_numaniso */) {
1258     NAMD_die("DID NOT FIND REQUIRED NUMANISO IN DRUDE PSF FILE");
1259   }
1260   else if (is_found_numaniso /* but not is_drude_file */) {
1261     NAMD_die("FOUND NUMANISO IN PSF FILE MISSING DRUDE DESIGNATION");
1262   }
1263
1264   if (is_found_ncrterm) {
1265     // We found crossterms section of PSF.
1266     // Read the number and then read the crossterms.
1267     sscanf(buffer, "%d", &numCrossterms);
1268     if (numCrossterms > 0) read_crossterms(psf_file, params);
1269   }
1270
1271   // Nothing else for us to read.
1272
1273   /*  Close the .psf file.  */
1274   Fclose(psf_file);
1275
1276   //  analyze the data and find the status of each atom
1277   numRealBonds = numBonds;
1278   build_atom_status();
1279   return;
1280 }
1281
1282 /************************************************************************/
1283 /*                  */
1284 /*        FUNCTION read_atoms      */
1285 /*                  */
1286 /*   INPUTS:                */
1287 /*  fd - file pointer to the .psf file        */
1288 /*  params - Parameters object to use for parameters    */
1289 /*                  */
1290 /*  this function reads in the Atoms section of the .psf file.      */
1291 /*   This section consists of numAtoms lines that are of the form:  */
1292 /*     <atom#> <mol> <seg#> <res> <atomname> <atomtype> <charge> <mass> */
1293 /*   Each line is read into the appropriate entry in the atoms array.   */
1294 /*   The parameters object is then used to determine the vdW constants  */
1295 /*   for this atom.              */
1296 /*                  */
1297 /************************************************************************/
1298
1299 void Molecule::read_atoms(FILE *fd, Parameters *params)
1300
1301 {
1302   char buffer[512];  // Buffer for reading from file
1303   int atom_number=0;  // Atom number 
1304   int last_atom_number=0; // Last atom number, used to assure
1305         // atoms are in order
1306   char segment_name[11]; // Segment name
1307   char residue_number[11]; // Residue number
1308   char residue_name[11];  // Residue name
1309   char atom_name[11];  // Atom name
1310   char atom_type[11];  // Atom type
1311   Real charge;    // Charge for the current atom
1312   Real mass;    // Mass for the current atom
1313   int read_count;    // Number of fields read by sscanf
1314
1315   /*  Allocate the atom arrays          */
1316   atoms     = new Atom[numAtoms];
1317   atomNames = new AtomNameInfo[numAtoms];
1318   if(simParams->genCompressedPsf) {
1319       atomSegResids = new AtomSegResInfo[numAtoms];
1320   }
1321
1322   // DRUDE: supplement Atom data
1323   if (is_drude_psf) {
1324     drudeConsts = new DrudeConst[numAtoms];
1325   }
1326   // DRUDE
1327
1328   hydrogenGroup.resize(0);
1329
1330   if (atoms == NULL || atomNames == NULL )
1331   {
1332     NAMD_die("memory allocation failed in Molecule::read_atoms");
1333   }
1334
1335   ResidueLookupElem *tmpResLookup = resLookup;
1336
1337   /*  Loop and read in numAtoms atom lines.      */
1338   while (atom_number < numAtoms)
1339   {
1340     // Standard PSF format has 8 columns:
1341     // ATOMNUM  SEGNAME  RESIDUE  RESNAME  ATOMNAME  ATOMTYPE  CHARGE  MASS
1342
1343     /*  Get the line from the file        */
1344     NAMD_read_line(fd, buffer);
1345
1346     /*  If its blank or a comment, skip it      */
1347     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') )
1348       continue;
1349
1350     /*  Parse up the line          */
1351     read_count=sscanf(buffer, "%d %s %s %s %s %s %f %f",
1352        &atom_number, segment_name, residue_number,
1353        residue_name, atom_name, atom_type, &charge, &mass);
1354
1355     /*  Check to make sure we found what we were expecting  */
1356     if (read_count != 8)
1357     {
1358       char err_msg[128];
1359
1360       sprintf(err_msg, "BAD ATOM LINE FORMAT IN PSF FILE IN ATOM LINE %d\nLINE=%s",
1361          last_atom_number+1, buffer);
1362       NAMD_die(err_msg);
1363     }
1364
1365     // DRUDE: read alpha and thole parameters from atom line
1366     if (is_drude_psf)
1367     {
1368       // Drude model PSF format has 11 columns, the 8 above plus 3 more:
1369       //   (unknown integer)  ALPHA  THOLE
1370       // These constants are used for the Thole interactions
1371       // (dipole interactions occurring between excluded non-bonded terms).
1372
1373       Real alpha, thole;
1374       read_count=sscanf(buffer,
1375 //          "%*d %*s %*s %*s %*s %*s %*f %*f %*d %*f %*f %f %f", &alpha, &thole);
1376                 // the two columns preceding alpha and thole will disappear
1377           "%*d %*s %*s %*s %*s %*s %*f %*f %*d %f %f", &alpha, &thole);
1378       if (read_count != 2)
1379       {
1380         char err_msg[128];
1381
1382         sprintf(err_msg, "BAD ATOM LINE FORMAT IN PSF FILE "
1383             "IN ATOM LINE %d\nLINE=%s", last_atom_number+1, buffer);
1384         NAMD_die(err_msg);
1385       }
1386       drudeConsts[atom_number-1].alpha = alpha;
1387       drudeConsts[atom_number-1].thole = thole;
1388     }
1389     // DRUDE
1390
1391     /*  Check if this is in XPLOR format  */
1392     int atom_type_num;
1393     if ( sscanf(atom_type, "%d", &atom_type_num) > 0 )
1394     {
1395       NAMD_die("Structure (psf) file is either in CHARMM format (with numbers for atoms types, the X-PLOR format using names is required) or the segment name field is empty.");
1396     }
1397
1398     /*  Make sure the atoms were in sequence    */
1399     if (atom_number != last_atom_number+1)
1400     {
1401       char err_msg[128];
1402
1403       sprintf(err_msg, "ATOM NUMBERS OUT OF ORDER AT ATOM #%d OF PSF FILE",
1404          last_atom_number+1);
1405       NAMD_die(err_msg);
1406     }
1407
1408     last_atom_number++;
1409
1410     /*  Dynamically allocate strings for atom name, atom    */
1411     /*  type, etc so that we only allocate as much space    */
1412     /*  for these strings as we really need      */
1413     int reslength = strlen(residue_name)+1;
1414     int namelength = strlen(atom_name)+1;
1415     int typelength = strlen(atom_type)+1;
1416
1417     atomNames[atom_number-1].resname = nameArena->getNewArray(reslength);
1418     atomNames[atom_number-1].atomname = nameArena->getNewArray(namelength);
1419     atomNames[atom_number-1].atomtype = nameArena->getNewArray(typelength);
1420   
1421     if (atomNames[atom_number-1].resname == NULL)
1422     {
1423       NAMD_die("memory allocation failed in Molecule::read_atoms");
1424     }
1425
1426     /*  Put the values from this atom into the atoms array  */
1427     strcpy(atomNames[atom_number-1].resname, residue_name);
1428     strcpy(atomNames[atom_number-1].atomname, atom_name);
1429     strcpy(atomNames[atom_number-1].atomtype, atom_type);
1430     atoms[atom_number-1].mass = mass;
1431     atoms[atom_number-1].charge = charge;
1432     atoms[atom_number-1].status = UnknownAtom;
1433
1434     /*  Add this atom to residue lookup table */
1435     if ( tmpResLookup ) tmpResLookup =
1436         tmpResLookup->append(segment_name, atoi(residue_number), atom_number-1);
1437
1438     if(atomSegResids) { //for compressing molecule information
1439         AtomSegResInfo *one = atomSegResids + (atom_number - 1);
1440         memcpy(one->segname, segment_name, strlen(segment_name)+1);
1441         one->resid = atoi(residue_number);
1442     }
1443
1444     /*  Determine the type of the atom (H or O) */
1445     if ( simParams->ignoreMass ) {
1446     } else if (atoms[atom_number-1].mass <= 0.05) {
1447       atoms[atom_number-1].status |= LonepairAtom;
1448     } else if (atoms[atom_number-1].mass < 1.0) {
1449       atoms[atom_number-1].status |= DrudeAtom;
1450     } else if (atoms[atom_number-1].mass <=3.5) {
1451       atoms[atom_number-1].status |= HydrogenAtom;
1452     } else if ((atomNames[atom_number-1].atomname[0] == 'O') && 
1453          (atoms[atom_number-1].mass >= 14.0) && 
1454          (atoms[atom_number-1].mass <= 18.0)) {
1455       atoms[atom_number-1].status |= OxygenAtom;
1456     }
1457
1458     /*  Look up the vdw constants for this atom    */
1459     params->assign_vdw_index(atomNames[atom_number-1].atomtype, 
1460        &(atoms[atom_number-1]));
1461         }
1462
1463   return;
1464 }
1465 /*      END OF FUNCTION read_atoms      */
1466
1467 /************************************************************************/
1468 /*                  */
1469 /*      FUNCTION read_bonds        */
1470 /*                  */
1471 /*  read_bonds reads in the bond section of the .psf file.  This    */
1472 /*  section contains a list of pairs of numbers where each pair is      */
1473 /*  represents two atoms that are bonded together.  Each atom pair is   */
1474 /*  read in.  Then that parameter object is queried to determine the    */
1475 /*  force constant and rest distance for the bond.      */
1476 /*                  */
1477 /************************************************************************/
1478
1479 void Molecule::read_bonds(FILE *fd, Parameters *params)
1480
1481 {
1482   int atom_nums[2];  // Atom indexes for the bonded atoms
1483   register int j;      // Loop counter
1484   int num_read=0;    // Number of bonds read so far
1485   int origNumBonds = numBonds;   // number of bonds in file header
1486
1487   /*  Allocate the array to hold the bonds      */
1488   bonds=new Bond[numBonds];
1489
1490   if (bonds == NULL)
1491   {
1492     NAMD_die("memory allocations failed in Molecule::read_bonds");
1493   }
1494
1495   /*  Loop through and read in all the bonds      */
1496   while (num_read < numBonds)
1497   {
1498     /*  Loop and read in the two atom indexes    */
1499     for (j=0; j<2; j++)
1500     {
1501       /*  Read the atom number from the file.         */
1502       /*  Subtract 1 to convert the index from the    */
1503       /*  1 to NumAtoms used in the file to the       */
1504       /*  0 to NumAtoms-1 that we need    */
1505       atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
1506
1507       /*  Check to make sure the index isn't too big  */
1508       if (atom_nums[j] >= numAtoms)
1509       {
1510         char err_msg[128];
1511
1512         sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1513         NAMD_die(err_msg);
1514       }
1515     }
1516
1517     /*  Assign the atom indexes to the array element  */
1518     Bond *b = &(bonds[num_read]);
1519     b->atom1=atom_nums[0];
1520     b->atom2=atom_nums[1];
1521
1522     /*  Query the parameter object for the constants for    */
1523     /*  this bond            */
1524     params->assign_bond_index(
1525         atomNames[atom_nums[0]].atomtype,
1526         atomNames[atom_nums[1]].atomtype,
1527         b);
1528
1529     /*  Make sure this isn't a fake bond meant for shake in x-plor.  */
1530     Real k, x0;
1531     params->get_bond_params(&k,&x0,b->bond_type);
1532     if (is_lonepairs_psf) {
1533       // need to retain Lonepair bonds for Drude
1534       if ( k == 0. && !is_lp(b->atom1) && !is_lp(b->atom2)) --numBonds;
1535       else ++num_read;
1536     }
1537     else {
1538       if ( k == 0. ) --numBonds;  // fake bond
1539       else ++num_read;  // real bond
1540     }
1541   }
1542
1543   /*  Tell user about our subterfuge  */
1544   if ( numBonds != origNumBonds ) {
1545     iout << iWARN << "Ignored " << origNumBonds - numBonds <<
1546             " bonds with zero force constants.\n" << endi;
1547     iout << iWARN <<
1548         "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
1549   }
1550
1551   return;
1552 }
1553 /*      END OF FUNCTION read_bonds      */
1554
1555 /************************************************************************/
1556 /*                  */
1557 /*      FUNCTION read_angles        */
1558 /*                  */
1559 /*   INPUTS:                */
1560 /*  fd - File descriptor for .psf file        */
1561 /*  params - Parameters object to query for parameters    */
1562 /*                  */
1563 /*  read_angles reads the angle parameters from the .psf file.      */
1564 /*   This section of the .psf file consists of a list of triplets of    */
1565 /*   atom indexes.  Each triplet represents three atoms connected via   */
1566 /*   an angle bond.  The parameter object is queried to obtain the      */
1567 /*   constants for each bond.            */
1568 /*                  */
1569 /************************************************************************/
1570
1571 void Molecule::read_angles(FILE *fd, Parameters *params)
1572
1573 {
1574   int atom_nums[3];  //  Atom numbers for the three atoms
1575   register int j;      //  Loop counter
1576   int num_read=0;    //  Number of angles read so far
1577   int origNumAngles = numAngles;  // Number of angles in file
1578   /*  Alloc the array of angles          */
1579   angles=new Angle[numAngles];
1580
1581   if (angles == NULL)
1582   {
1583     NAMD_die("memory allocation failed in Molecule::read_angles");
1584   }
1585
1586   /*  Loop through and read all the angles      */
1587   while (num_read < numAngles)
1588   {
1589     /*  Loop through the 3 atom indexes in the current angle*/
1590     for (j=0; j<3; j++)
1591     {
1592       /*  Read the atom number from the file.         */
1593       /*  Subtract 1 to convert the index from the    */
1594       /*  1 to NumAtoms used in the file to the       */
1595       /*  0 to NumAtoms-1 that we need    */
1596       atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
1597
1598       /*  Check to make sure the atom index doesn't   */
1599       /*  exceed the Number of Atoms      */
1600       if (atom_nums[j] >= numAtoms)
1601       {
1602         char err_msg[128];
1603
1604         sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1605         NAMD_die(err_msg);
1606       }
1607     }
1608
1609     /*  Assign the three atom indices      */
1610     angles[num_read].atom1=atom_nums[0];
1611     angles[num_read].atom2=atom_nums[1];
1612     angles[num_read].atom3=atom_nums[2];
1613
1614     /*  Get the constant values for this bond from the  */
1615     /*  parameter object          */
1616     params->assign_angle_index(
1617       atomNames[atom_nums[0]].atomtype,
1618       atomNames[atom_nums[1]].atomtype,
1619       atomNames[atom_nums[2]].atomtype,
1620       &(angles[num_read]), simParams->alchOn ? -1 : 0);
1621     if ( angles[num_read].angle_type == -1 ) {
1622       iout << iWARN << "ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
1623            << endi;
1624     }
1625
1626     /*  Make sure this isn't a fake angle meant for shake in x-plor.  */
1627     Real k, t0, k_ub, r_ub;
1628     if ( angles[num_read].angle_type == -1 ) { k = -1.;  k_ub = -1.; } else
1629     params->get_angle_params(&k,&t0,&k_ub,&r_ub,angles[num_read].angle_type);
1630     if ( k == 0. && k_ub == 0. ) --numAngles;  // fake angle
1631     else ++num_read;  // real angle
1632   }
1633
1634   /*  Tell user about our subterfuge  */
1635   if ( numAngles != origNumAngles ) {
1636     iout << iWARN << "Ignored " << origNumAngles - numAngles <<
1637             " angles with zero force constants.\n" << endi;
1638   }
1639
1640   return;
1641 }
1642 /*      END OF FUNCTION read_angles      */
1643
1644 /************************************************************************/
1645 /*                  */
1646 /*        FUNCTION read_dihedrals      */
1647 /*                  */
1648 /*   INPUTS:                */
1649 /*  fd - file descriptor for the .psf file        */
1650 /*  params - pointer to parameter object        */
1651 /*                  */
1652 /*  read_dihedreals reads the dihedral section of the .psf file.    */
1653 /*   This section of the file contains a list of quartets of atom       */
1654 /*   numbers.  Each quartet represents a group of atoms that form a     */
1655 /*   dihedral bond.              */
1656 /*                  */
1657 /************************************************************************/
1658
1659 void Molecule::read_dihedrals(FILE *fd, Parameters *params)
1660 {
1661   int atom_nums[4];  // The 4 atom indexes
1662   int last_atom_nums[4];  // Atom numbers from previous bond
1663   register int j;      // loop counter
1664   int num_read=0;    // number of dihedrals read so far
1665   int multiplicity=1;  // multiplicity of the current bond
1666   Bool duplicate_bond;  // Is this a duplicate of the last bond
1667   int num_unique=0;   // Number of unique dihedral bonds
1668
1669   //  Initialize the array used to check for duplicate dihedrals
1670   for (j=0; j<4; j++)
1671     last_atom_nums[j] = -1;
1672
1673   /*  Allocate an array to hold the Dihedrals      */
1674   dihedrals = new Dihedral[numDihedrals];
1675
1676   if (dihedrals == NULL)
1677   {
1678     NAMD_die("memory allocation failed in Molecule::read_dihedrals");
1679   }
1680
1681   /*  Loop through and read all the dihedrals      */
1682   while (num_read < numDihedrals)
1683   {
1684     duplicate_bond = TRUE;
1685
1686     /*  Loop through and read the 4 indexes for this bond   */
1687     for (j=0; j<4; j++)
1688     {
1689       /*  Read the atom number from the file.         */
1690       /*  Subtract 1 to convert the index from the    */
1691       /*  1 to NumAtoms used in the file to the       */
1692       /*  0 to NumAtoms-1 that we need    */
1693       atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
1694
1695       /*  Check for an atom index that is too large  */
1696       if (atom_nums[j] >= numAtoms)
1697       {
1698         char err_msg[128];
1699
1700         sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1701         NAMD_die(err_msg);
1702       }
1703
1704       //  Check to see if this atom matches the last bond
1705       if (atom_nums[j] != last_atom_nums[j])
1706       {
1707         duplicate_bond = FALSE;
1708       }
1709
1710       last_atom_nums[j] = atom_nums[j];
1711     }
1712
1713     //  Check to see if this is really a new bond or just
1714     //  a repeat of the last one
1715     if (duplicate_bond)
1716     {
1717       //  This is a duplicate, so increase the multiplicity
1718       multiplicity++;
1719
1720       if (multiplicity == 2)
1721       {
1722         numMultipleDihedrals++;
1723       }
1724     }
1725     else
1726     {
1727       multiplicity=1;
1728       num_unique++;
1729     }
1730
1731     /*  Assign the atom indexes        */
1732     dihedrals[num_unique-1].atom1=atom_nums[0];
1733     dihedrals[num_unique-1].atom2=atom_nums[1];
1734     dihedrals[num_unique-1].atom3=atom_nums[2];
1735     dihedrals[num_unique-1].atom4=atom_nums[3];
1736
1737     /*  Get the constants for this dihedral bond    */
1738     params->assign_dihedral_index(
1739        atomNames[atom_nums[0]].atomtype,
1740        atomNames[atom_nums[1]].atomtype,
1741        atomNames[atom_nums[2]].atomtype,
1742        atomNames[atom_nums[3]].atomtype,
1743        &(dihedrals[num_unique-1]),
1744        multiplicity, simParams->alchOn ? -1 : 0);
1745     if ( dihedrals[num_unique-1].dihedral_type == -1 ) {
1746       iout << iWARN << "ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
1747            << endi;
1748     }
1749
1750     num_read++;
1751   }
1752
1753   numDihedrals = num_unique;
1754
1755   return;
1756 }
1757 /*      END OF FUNCTION read_dihedral      */
1758
1759 /************************************************************************/
1760 /*                  */
1761 /*        FUNCTION read_impropers      */
1762 /*                  */
1763 /*   INPUTS:                */
1764 /*  fd - file descriptor for .psf file        */
1765 /*  params - parameter object          */
1766 /*                  */
1767 /*  read_impropers reads the improper section of the .psf file.  */
1768 /*   This section is identical to the dihedral section in that it is    */
1769 /*   made up of a list of quartets of atom indexes that define the      */
1770 /*   atoms that are bonded together.          */
1771 /*                  */
1772 /************************************************************************/
1773
1774 void Molecule::read_impropers(FILE *fd, Parameters *params)
1775 {
1776   int atom_nums[4];  //  Atom indexes for the 4 atoms
1777   int last_atom_nums[4];  //  Atom indexes from previous bond
1778   register int j;      //  Loop counter
1779   int num_read=0;    //  Number of impropers read so far
1780   int multiplicity=1;  // multiplicity of the current bond
1781   Bool duplicate_bond;  // Is this a duplicate of the last bond
1782   int num_unique=0;   // Number of unique dihedral bonds
1783
1784   //  Initialize the array used to look for duplicate improper
1785   //  entries.  Set them all to -1 so we know nothing will match
1786   for (j=0; j<4; j++)
1787     last_atom_nums[j] = -1;
1788
1789   /*  Allocate the array to hold the impropers      */
1790   impropers=new Improper[numImpropers];
1791
1792   if (impropers == NULL)
1793   {
1794     NAMD_die("memory allocation failed in Molecule::read_impropers");
1795   }
1796
1797   /*  Loop through and read all the impropers      */
1798   while (num_read < numImpropers)
1799   {
1800     duplicate_bond = TRUE;
1801
1802     /*  Loop through the 4 indexes for this improper  */
1803     for (j=0; j<4; j++)
1804     {
1805       /*  Read the atom number from the file.         */
1806       /*  Subtract 1 to convert the index from the    */
1807       /*  1 to NumAtoms used in the file to the       */
1808       /*  0 to NumAtoms-1 that we need    */
1809       atom_nums[j]=NAMD_read_int(fd, "IMPROPERS")-1;
1810
1811       /*  Check to make sure the index isn't too big  */
1812       if (atom_nums[j] >= numAtoms)
1813       {
1814         char err_msg[128];
1815
1816         sprintf(err_msg, "IMPROPERS INDEX %d GREATER THAN NATOM %d IN IMPROPERS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1817         NAMD_die(err_msg);
1818       }
1819
1820       if (atom_nums[j] != last_atom_nums[j])
1821       {
1822         duplicate_bond = FALSE;
1823       }
1824
1825       last_atom_nums[j] = atom_nums[j];
1826     }
1827
1828     //  Check to see if this is a duplicate improper
1829     if (duplicate_bond)
1830     {
1831       //  This is a duplicate improper.  So we don't
1832       //  really count this entry, we just update
1833       //  the parameters object
1834       multiplicity++;
1835
1836       if (multiplicity == 2)
1837       {
1838         //  Count the number of multiples.
1839         numMultipleImpropers++;
1840       }
1841     }
1842     else
1843     {
1844       //  Not a duplicate
1845       multiplicity = 1;
1846       num_unique++;
1847     }
1848
1849     /*  Assign the atom indexes        */
1850     impropers[num_unique-1].atom1=atom_nums[0];
1851     impropers[num_unique-1].atom2=atom_nums[1];
1852     impropers[num_unique-1].atom3=atom_nums[2];
1853     impropers[num_unique-1].atom4=atom_nums[3];
1854
1855     /*  Look up the constants for this bond      */
1856     params->assign_improper_index(
1857        atomNames[atom_nums[0]].atomtype,
1858        atomNames[atom_nums[1]].atomtype,
1859        atomNames[atom_nums[2]].atomtype,
1860        atomNames[atom_nums[3]].atomtype,
1861        &(impropers[num_unique-1]),
1862        multiplicity);
1863
1864     num_read++;
1865   }
1866
1867   //  Now reset the numImpropers value to the number of UNIQUE
1868   //  impropers.  Sure, we waste a few entries in the improper_array
1869   //  on the master node, but it is very little space . . .
1870   numImpropers = num_unique;
1871
1872   return;
1873 }
1874 /*      END OF FUNCTION read_impropers      */
1875
1876 /************************************************************************/
1877 /*                  */
1878 /*        FUNCTION read_crossterms      */
1879 /*                  */
1880 /*   INPUTS:                */
1881 /*  fd - file descriptor for .psf file        */
1882 /*  params - parameter object          */
1883 /*                  */
1884 /*   This section is identical to the dihedral section in that it is    */
1885 /*   made up of a list of quartets of atom indexes that define the      */
1886 /*   atoms that are bonded together.          */
1887 /*                  */
1888 /************************************************************************/
1889
1890 void Molecule::read_crossterms(FILE *fd, Parameters *params)
1891
1892 {
1893   int atom_nums[8];  //  Atom indexes for the 4 atoms
1894   int last_atom_nums[8];  //  Atom indexes from previous bond
1895   register int j;      //  Loop counter
1896   int num_read=0;    //  Number of items read so far
1897   Bool duplicate_bond;  // Is this a duplicate of the last bond
1898
1899   //  Initialize the array used to look for duplicate crossterm
1900   //  entries.  Set them all to -1 so we know nothing will match
1901   for (j=0; j<8; j++)
1902     last_atom_nums[j] = -1;
1903
1904   /*  Allocate the array to hold the cross-terms */
1905   crossterms=new Crossterm[numCrossterms];
1906
1907   if (crossterms == NULL)
1908   {
1909     NAMD_die("memory allocation failed in Molecule::read_crossterms");
1910   }
1911
1912   /*  Loop through and read all the cross-terms      */
1913   while (num_read < numCrossterms)
1914   {
1915     duplicate_bond = TRUE;
1916
1917     /*  Loop through the 8 indexes for this cross-term */
1918     for (j=0; j<8; j++)
1919     {
1920       /*  Read the atom number from the file.         */
1921       /*  Subtract 1 to convert the index from the    */
1922       /*  1 to NumAtoms used in the file to the       */
1923       /*  0 to NumAtoms-1 that we need    */
1924       atom_nums[j]=NAMD_read_int(fd, "CROSS-TERMS")-1;
1925
1926       /*  Check to make sure the index isn't too big  */
1927       if (atom_nums[j] >= numAtoms)
1928       {
1929         char err_msg[128];
1930
1931         sprintf(err_msg, "CROSS-TERM INDEX %d GREATER THAN NATOM %d IN CROSS-TERMS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1932         NAMD_die(err_msg);
1933       }
1934
1935       if (atom_nums[j] != last_atom_nums[j])
1936       {
1937         duplicate_bond = FALSE;
1938       }
1939
1940       last_atom_nums[j] = atom_nums[j];
1941     }
1942
1943     //  Check to see if this is a duplicate term
1944     if (duplicate_bond)
1945     {
1946       iout << iWARN << "Duplicate cross-term detected.\n" << endi;
1947     }
1948
1949     /*  Assign the atom indexes        */
1950     crossterms[num_read].atom1=atom_nums[0];
1951     crossterms[num_read].atom2=atom_nums[1];
1952     crossterms[num_read].atom3=atom_nums[2];
1953     crossterms[num_read].atom4=atom_nums[3];
1954     crossterms[num_read].atom5=atom_nums[4];
1955     crossterms[num_read].atom6=atom_nums[5];
1956     crossterms[num_read].atom7=atom_nums[6];
1957     crossterms[num_read].atom8=atom_nums[7];
1958
1959     /*  Look up the constants for this bond      */
1960     params->assign_crossterm_index(
1961        atomNames[atom_nums[0]].atomtype,
1962        atomNames[atom_nums[1]].atomtype,
1963        atomNames[atom_nums[2]].atomtype,
1964        atomNames[atom_nums[3]].atomtype,
1965        atomNames[atom_nums[4]].atomtype,
1966        atomNames[atom_nums[5]].atomtype,
1967        atomNames[atom_nums[6]].atomtype,
1968        atomNames[atom_nums[7]].atomtype,
1969        &(crossterms[num_read]));
1970
1971     if(!duplicate_bond) num_read++;
1972   }
1973
1974   numCrossterms = num_read;
1975
1976   return;
1977 }
1978 /*      END OF FUNCTION read_impropers      */
1979
1980 /************************************************************************/
1981 /*                  */
1982 /*      FUNCTION read_donors        */
1983 /*                  */
1984 /*  read_donors reads in the bond section of the .psf file.  This   */
1985 /*  section contains a list of pairs of numbers where each pair is      */
1986 /*  represents two atoms that are part of an H-bond.  Each atom pair is */
1987 /*  read in.                                                            */
1988 /*                  */
1989 /*  Donor atoms are the heavy atoms to which hydrogens are bonded.      */
1990 /*  There will always be a donor atom for each donor pair.  However,    */
1991 /*  for a united-atom model there may not be an explicit hydrogen       */
1992 /*  present, in which case the second atom index in the pair will be    */
1993 /*  given as 0 in the PSF (and stored as -1 in this program's internal  */
1994 /*  storage).                                                           */
1995 /************************************************************************/
1996
1997 void Molecule::read_donors(FILE *fd)
1998 {
1999   int d[2];               // temporary storage of donor atom index
2000   register int j;      // Loop counter
2001   int num_read=0;    // Number of bonds read so far
2002   int num_no_hydr=0;      // Number of bonds with no hydrogen given
2003
2004   /*  Allocate the array to hold the bonds      */
2005   donors=new Bond[numDonors];
2006
2007   if (donors == NULL)
2008   {
2009     NAMD_die("memory allocations failed in Molecule::read_donors");
2010   }
2011
2012   /*  Loop through and read in all the donors      */
2013   while (num_read < numDonors)
2014   {
2015     /*  Loop and read in the two atom indexes    */
2016     for (j=0; j<2; j++)
2017     {
2018       /*  Read the atom number from the file.         */
2019       /*  Subtract 1 to convert the index from the    */
2020       /*  1 to NumAtoms used in the file to the       */
2021       /*  0 to NumAtoms-1 that we need    */
2022       d[j]=NAMD_read_int(fd, "DONORS")-1;
2023
2024       /*  Check to make sure the index isn't too big  */
2025       if (d[j] >= numAtoms)
2026       {
2027         char err_msg[128];
2028
2029         sprintf(err_msg,
2030     "DONOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE",
2031           d[j]+1, numAtoms, num_read+1);
2032         NAMD_die(err_msg);
2033       }
2034
2035       /*  Check if there is a hydrogen given */
2036       if (d[j] < 0)
2037                           num_no_hydr++;
2038     }
2039
2040     /*  Assign the atom indexes to the array element  */
2041     Bond *b = &(donors[num_read]);
2042     b->atom1=d[0];
2043     b->atom2=d[1];
2044
2045     num_read++;
2046   }
2047
2048   return;
2049 }
2050 /*      END OF FUNCTION read_donors      */
2051
2052
2053 /************************************************************************/
2054 /*                  */
2055 /*      FUNCTION read_acceptors        */
2056 /*                  */
2057 /*  read_acceptors reads in the bond section of the .psf file.      */
2058 /*  This section contains a list of pairs of numbers where each pair is */
2059 /*  represents two atoms that are part of an H-bond.  Each atom pair is */
2060 /*  read in.                                                            */
2061 /*                  */
2062 /*  Acceptor atoms are the heavy atoms to which hydrogens directly      */
2063 /*  orient in a hydrogen bond interaction.  There will always be an     */
2064 /*  acceptor atom for each acceptor pair.  The antecedent atom, to      */
2065 /*  which the acceptor is bound, may not be given in the structure,     */
2066 /*  however, in which case the second atom index in the pair will be    */
2067 /*  given as 0 in the PSF (and stored as -1 in this program's internal  */
2068 /*  storage).                                                           */
2069 /************************************************************************/
2070
2071 void Molecule::read_acceptors(FILE *fd)
2072 {
2073   int d[2];               // temporary storage of atom index
2074   register int j;      // Loop counter
2075   int num_read=0;    // Number of bonds read so far
2076         int num_no_ante=0;      // number of pairs with no antecedent
2077
2078   /*  Allocate the array to hold the bonds      */
2079   acceptors=new Bond[numAcceptors];
2080
2081   if (acceptors == NULL)
2082   {
2083     NAMD_die("memory allocations failed in Molecule::read_acceptors");
2084   }
2085
2086   /*  Loop through and read in all the acceptors      */
2087   while (num_read < numAcceptors)
2088   {
2089     /*  Loop and read in the two atom indexes    */
2090     for (j=0; j<2; j++)
2091     {
2092       /*  Read the atom number from the file.         */
2093       /*  Subtract 1 to convert the index from the    */
2094       /*  1 to NumAtoms used in the file to the       */
2095       /*  0 to NumAtoms-1 that we need    */
2096       d[j]=NAMD_read_int(fd, "ACCEPTORS")-1;
2097
2098       /*  Check to make sure the index isn't too big  */
2099       if (d[j] >= numAtoms)
2100       {
2101         char err_msg[128];
2102
2103         sprintf(err_msg, "ACCEPTOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE", d[j]+1, numAtoms, num_read+1);
2104         NAMD_die(err_msg);
2105       }
2106
2107       /*  Check if there is an antecedent given */
2108       if (d[j] < 0)
2109                           num_no_ante++;
2110     }
2111
2112     /*  Assign the atom indexes to the array element  */
2113     Bond *b = &(acceptors[num_read]);
2114     b->atom1=d[0];
2115     b->atom2=d[1];
2116
2117     num_read++;
2118   }
2119
2120   return;
2121 }
2122 /*      END OF FUNCTION read_acceptors      */
2123
2124
2125 /************************************************************************/
2126 /*                  */
2127 /*      FUNCTION read_exclusions      */
2128 /*                  */
2129 /*   INPUTS:                */
2130 /*  fd - file descriptor for .psf file        */
2131 /*                  */
2132 /*  read_exclusions reads in the explicit non-bonded exclusions     */
2133 /*  from the .psf file.  This section is a little funky, so hang on.    */
2134 /*  Ok, first there is a list of atom indexes that is NumExclusions     */
2135 /*  long.  These are in some sense the atoms that will be exlcuded.     */
2136 /*  Following this list is a list of NumAtoms length that is a list     */
2137 /*  of indexes into the list of excluded atoms.  So an example.  Suppose*/
2138 /*  we have a 5 atom simulation with 3 explicit exclusions.  The .psf   */
2139 /*  file could look like:            */
2140 /*                  */
2141 /*  3!NNB                */
2142 /*  3 4 5                */
2143 /*  0 1 3 3 3              */
2144 /*                  */
2145 /*  This would mean that atom 1 has no explicit exclusions.  Atom 2     */
2146 /*  has an explicit exclusion with atom 3.  Atom 3 has an explicit      */
2147 /*  exclusion with atoms 4 AND 5.  And atoms 4 and 5 have no explicit   */
2148 /*  exclusions.  Got it!?!  I'm not sure who dreamed this up . . .      */
2149 /*                  */
2150 /************************************************************************/
2151
2152 void Molecule::read_exclusions(FILE *fd)
2153
2154 {
2155   int *exclusion_atoms;  //  Array of indexes of excluded atoms
2156   register int num_read=0;    //  Number fo exclusions read in
2157   int current_index;  //  Current index value
2158   int last_index;    //  the previous index value
2159   register int insert_index=0;  //  index of where we are in exlcusions array
2160
2161   /*  Allocate the array of exclusion structures and the array of */
2162   /*  exlcuded atom indexes          */
2163   exclusions      = new Exclusion[numExclusions];
2164   exclusion_atoms = new int[numExclusions];
2165
2166   if ( (exclusions == NULL) || (exclusion_atoms == NULL) )
2167   {
2168     NAMD_die("memory allocation failed in Molecule::read_exclusions");
2169   }
2170
2171   /*  First, read in the excluded atoms list      */
2172   for (num_read=0; num_read<numExclusions; num_read++)
2173   {
2174     /*  Read the atom number from the file. Subtract 1 to   */
2175     /*  convert the index from the 1 to NumAtoms used in the*/
2176     /*  file to the  0 to NumAtoms-1 that we need    */
2177     exclusion_atoms[num_read]=NAMD_read_int(fd, "IMPROPERS")-1;
2178
2179     /*  Check for an illegal index        */
2180     if (exclusion_atoms[num_read] >= numAtoms)
2181     {
2182       char err_msg[128];
2183
2184       sprintf(err_msg, "EXCLUSION INDEX %d GREATER THAN NATOM %d IN EXCLUSION # %d IN PSF FILE", exclusion_atoms[num_read]+1, numAtoms, num_read+1);
2185       NAMD_die(err_msg);
2186     }
2187   }
2188
2189   /*  Now, go through and read the list of NumAtoms pointers into */
2190   /*  the array that we just read in        */
2191   last_index=0;
2192
2193   for (num_read=0; num_read<numAtoms; num_read++)
2194   {
2195     /*  Read in the current index value      */
2196     current_index=NAMD_read_int(fd, "EXCLUSIONS");
2197
2198     /*  Check for an illegal pointer      */
2199     if (current_index>numExclusions)
2200     {
2201       char err_msg[128];
2202
2203       sprintf(err_msg, "EXCLUSION INDEX %d LARGER THAN NUMBER OF EXLCUSIONS %d IN PSF FILE, EXCLUSION #%d\n", 
2204          current_index+1, numExclusions, num_read);
2205       NAMD_die(err_msg);
2206     }
2207
2208     /*  Check to see if it matches the last index.  If so   */
2209     /*  than this atom has no exclusions.  If not, then     */
2210     /*  we have to build some exclusions      */
2211     if (current_index != last_index)
2212     {
2213       /*  This atom has some exclusions.  Loop from   */
2214       /*  the last_index to the current index.  This  */
2215       /*  will include how ever many exclusions this  */
2216       /*  atom has          */
2217       for (insert_index=last_index; 
2218            insert_index<current_index; insert_index++)
2219       {
2220         /*  Assign the two atoms involved.      */
2221         /*  The first one is our position in    */
2222         /*  the list, the second is based on    */
2223         /*  the pointer into the index list     */
2224         int a1 = num_read;
2225         int a2 = exclusion_atoms[insert_index];
2226         if ( a1 < a2 ) {
2227           exclusions[insert_index].atom1 = a1;
2228           exclusions[insert_index].atom2 = a2;
2229         } else if ( a2 < a1 ) {
2230           exclusions[insert_index].atom1 = a2;
2231           exclusions[insert_index].atom2 = a1;
2232         } else {
2233           char err_msg[128];
2234           sprintf(err_msg, "ATOM %d EXCLUDED FROM ITSELF IN PSF FILE\n", a1+1);
2235           NAMD_die(err_msg);
2236         }
2237       }
2238
2239       last_index=current_index;
2240     }
2241   }
2242
2243   /*  Free our temporary list of indexes        */
2244   delete [] exclusion_atoms;
2245
2246   return;
2247 }
2248 /*      END OF FUNCTION read_exclusions      */
2249
2250 /************************************************************************/
2251 /*      FUNCTION read_exclusions      */
2252 /*                  */
2253 /*   INPUTS:                */
2254 /*   int* atom_i - array of atom i indices    */
2255 /*   int* atom_j - array of atom j indices    */
2256 /*   int num_exclusion - length of array           */
2257 /*                  */
2258 /* JLai August 16th, 2012 */
2259 /************************************************************************/
2260 void Molecule::read_exclusions(int* atom_i, int* atom_j, int num_exclusion)
2261 {
2262     /*  Allocate the array of exclusion structures and the array of */
2263   /*  exlcuded atom indexes          */
2264   exclusions       = new Exclusion[num_exclusion];
2265   int loop_counter = 0;  
2266   int a=0;
2267   int b=0;
2268
2269   if ( (exclusions == NULL) )
2270   {
2271     NAMD_die("memory allocation failed in Molecule::read_exclusions");
2272   }
2273
2274   /* The following code only guarantees that exclusion.atom1 is < exclusion.atom2 */
2275   for (loop_counter = 0; loop_counter < num_exclusion; loop_counter++) {
2276         
2277         if ( (atom_i == NULL) || (atom_j == NULL) ) {
2278           NAMD_die("null pointer expection in Molecule::read_exclusions");
2279         }
2280
2281         a = atom_i[loop_counter];
2282         b = atom_j[loop_counter];
2283         if(a < b) {
2284                 exclusions[loop_counter].atom1 = a;
2285                 exclusions[loop_counter].atom2 = b;
2286         } else {
2287                 exclusions[loop_counter].atom1 = b;
2288                 exclusions[loop_counter].atom2 = a;
2289         }
2290         exclusionSet.add(Exclusion(exclusions[loop_counter].atom1,exclusions[loop_counter].atom2));
2291   }
2292
2293   if ( ! CkMyPe() ) {
2294     iout << iINFO << "ADDED " << num_exclusion << " EXPLICIT EXCLUSIONS: THIS VALUE WILL *NOT* BE ADDED TO THE STRUCTURE SUMMARY\n" << endi;
2295   }
2296
2297    return;
2298 }
2299 /*      END OF FUNCTION read_exclusions      */
2300
2301 /************************************************************************
2302  *
2303  *        FUNCTION read_lphosts
2304  *
2305  *   INPUTS:
2306  *  fd - file pointer to the .psf file
2307  *
2308  *  This function reads in the lone pair host section of the .psf file.
2309  *  Each lone pair is supported by 2 or 3 host atoms.
2310  *
2311  *  All lonepair specifications in a PSF are expected to have three
2312  *  associated values.  However, the meaning of these values depends
2313  *  on the lonepair type.  Nonetheless, for simplicity with the old
2314  *  code, the struct values are still called "distance", "angle",
2315  *  and "dihedral", even when this is not what the value signifies.
2316  *
2317  ************************************************************************/
2318 void Molecule::read_lphosts(FILE *fd)
2319 {
2320   char buffer[512];  // Buffer for reading from file
2321   char weight[8];    // Weighting type identified by string --
2322                      // unused/unsupported outside of CHARMM RTF
2323                      // so we read it, but ignore it.
2324   int numhosts;  // Refers to the number of atoms that support the
2325                  // given lone pair, either 2 or 3.
2326   int index;  // 1-based index into the stream of numbers (8 per line)
2327               // that indicates the start of each LP host record.
2328   int next_index = 1;  // Forecast next index value as an error check.
2329   int i, read_count;
2330   Real value1, value2, value3; 
2331
2332   // called only if numLphosts > 0
2333   lphosts = new Lphost[numLphosts];
2334   if (lphosts == NULL) {
2335     NAMD_die("memory allocation failed in Molecule::read_lphosts");
2336   }
2337   for (i = 0;  i < numLphosts;  i++) {
2338     NAMD_read_line(fd, buffer);
2339     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
2340     read_count=sscanf(buffer, "%d %d %6s %f %f %f",
2341         &numhosts, &index, weight, &value1, &value2, &value3);
2342     // The "weight" specification in PSF remains hardcoded as "F"
2343     // (for false) inside NAMD.  Again, no extant force field uses
2344     // lonepairs with weighted placement, so this can really only be
2345     // specified inside an RTF, which NAMD never reads anyway.
2346     if (read_count != 6 || index != next_index ||
2347         strcmp(weight,"F") != 0 || numhosts < 2 || 3 < numhosts) {
2348       char err_msg[128];
2349       sprintf(err_msg, "BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n"
2350           "LINE=%s\n", i+1, buffer);
2351       NAMD_die(err_msg);
2352     }
2353     lphosts[i].numhosts = numhosts; // currently must be 2 or 3
2354     next_index += numhosts + 1;  // add 1 to account for LP index
2355     if (numhosts == 2) {
2356       lphosts[i].distance = value1;
2357       lphosts[i].angle = value2;
2358       lphosts[i].dihedral = 0.0;  // ignore value3
2359     }
2360     else {  // numhosts == 3
2361       lphosts[i].distance = value1;
2362       lphosts[i].angle = value2 * (M_PI/180);     // convert to radians
2363       lphosts[i].dihedral = value3 * (M_PI/180);  // convert to radians
2364     }
2365   }
2366   for (i = 0;  i < numLphosts;  i++) {
2367     // Subtract 1 to get 0-based atom index
2368     lphosts[i].atom1 = NAMD_read_int(fd, "LPHOSTS")-1;
2369     lphosts[i].atom2 = NAMD_read_int(fd, "LPHOSTS")-1;
2370     lphosts[i].atom3 = NAMD_read_int(fd, "LPHOSTS")-1;
2371     // For numhosts==2, set unused atom4 to atom1
2372     lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
2373         NAMD_read_int(fd, "LPHOSTS")-1 : lphosts[i].atom1 );
2374   }
2375 }
2376 /*      END OF FUNCTION read_lphosts    */
2377
2378 /************************************************************************/
2379 /*                  */
2380 /*        FUNCTION read_anisos     */
2381 /*                  */
2382 /*   INPUTS:                */
2383 /*  fd - file pointer to the .psf file        */
2384 /*                  */
2385 /*  this function reads in the anisotropic terms section of .psf file. */
2386 /*                  */
2387 void Molecule::read_anisos(FILE *fd)
2388 {
2389   char buffer[512];  // Buffer for reading from file
2390   int numhosts, index, i, read_count;
2391   Real k11, k22, k33;
2392
2393   anisos = new Aniso[numAnisos];
2394   if (anisos == NULL)
2395   {
2396     NAMD_die("memory allocation failed in Molecule::read_anisos");
2397   }
2398   for (i = 0;  i < numAnisos;  i++)
2399   {
2400     NAMD_read_line(fd, buffer);
2401     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
2402     read_count=sscanf(buffer, "%f %f %f", &k11, &k22, &k33);
2403     if (read_count != 3)
2404     {
2405       char err_msg[128];
2406       sprintf(err_msg, "BAD FORMAT FOR ANISO LINE %d IN PSF FILE LINE\n"
2407           "LINE=%s\n", i+1, buffer);
2408       NAMD_die(err_msg);
2409     }
2410     anisos[i].k11 = k11;
2411     anisos[i].k22 = k22;
2412     anisos[i].k33 = k33;
2413   }
2414   for (i = 0;  i < numAnisos;  i++) {
2415     anisos[i].atom1 = NAMD_read_int(fd, "ANISOS")-1;
2416     anisos[i].atom2 = NAMD_read_int(fd, "ANISOS")-1;
2417     anisos[i].atom3 = NAMD_read_int(fd, "ANISOS")-1;
2418     anisos[i].atom4 = NAMD_read_int(fd, "ANISOS")-1;
2419   }
2420 }
2421 /*      END OF FUNCTION read_anisos     */
2422
2423 //LCPO
2424 inline int getLCPOTypeAmber(char atomType[11], int numBonds) {
2425
2426   //Hydrogen
2427   if (atomType[0] == 'H' || atomType[0] == 'h') {
2428     return 0;
2429
2430   //Carbon
2431   } else if (atomType[0] == 'C' || atomType[0] == 'c') {
2432     if (//Sp3 Carbon
2433       //atomType[1] == 'T')// ||
2434       strcmp(atomType, "CT" )==0 )
2435       //strcmp(atomType, "CP1" )==0 ||
2436       //strcmp(atomType, "CP2" )==0 ||
2437       //strcmp(atomType, "CP3" )==0 ||
2438       //strcmp(atomType, "CS"  )==0 )
2439       {
2440       if (numBonds == 1)
2441         return 1;
2442       else if (numBonds == 2)
2443         return 2;
2444       else if (numBonds == 3)
2445         return 3;
2446       else if (numBonds == 4)
2447         return 4;
2448       else
2449         return 1;
2450
2451     } else {//Sp2 or other
2452       if (numBonds == 2)
2453         return 5;
2454       else if (numBonds == 3)
2455         return 6;
2456       else
2457         return 1;
2458     }
2459
2460   //Nitrogen
2461   } else if (atomType[0] == 'N' || atomType[0] == 'n') {
2462     if ( strcmp(atomType, "N3"  ) == 0 ) { //Sp3 Nitrogen
2463       if (numBonds == 1)
2464         return 11;
2465       else if (numBonds == 2)
2466         return 12;
2467       else if (numBonds == 3)
2468         return 13;
2469       else
2470         return 11;
2471
2472     } else {//SP2 Nitrogen
2473       if (numBonds == 1)
2474         return 14;
2475       else if (numBonds == 2)
2476         return 15;
2477       else if (numBonds == 3)
2478         return 16;
2479       else
2480         return 11;
2481     }
2482
2483   //Oxygen
2484   } else if (atomType[0] == 'O' || atomType[0] == 'o') {
2485
2486     if ( strcmp(atomType, "O" )==0) {//Sp2 Oxygen
2487       return 9;
2488     } else if (strcmp(atomType, "O2" )==0) {//Carboxylate Oxygen
2489       return 10;
2490     } else { // Sp3 Oxygen
2491       if (numBonds == 1)
2492         return 7;
2493       else if (numBonds == 2)
2494         return 8;
2495       else
2496         return 7;
2497     }
2498
2499   //Sulfur
2500   } else if (atomType[0] == 'S' || atomType[0] == 's') {
2501     if ( strcmp(atomType, "SH" )==0) { //Sulfur 1 neighbor
2502       return 17;
2503     } else {
2504       return 18;
2505     }
2506
2507   //Phosphorus
2508   } else if (atomType[0] == 'P' || atomType[0] == 'p') {
2509       if (numBonds == 3)
2510         return 19;
2511       else if (numBonds == 4)
2512         return 20;
2513       else
2514         return 19;
2515   } else if (atomType[0] == 'Z') { // ? just to agree with Amber mdread.f
2516     return 0;
2517   } else  if ( strcmp(atomType, "MG" )==0) { //Mg
2518     return 22;
2519   } else { // unknown atom type
2520     return 5;
2521   }
2522   return 5;
2523 } // getLCPOTypeAmber
2524
2525 inline int getLCPOTypeCharmm(char atomType[11], int numBonds) {
2526
2527   //Hydrogen
2528   if (atomType[0] == 'H') {
2529     return 0;
2530
2531   //Carbon
2532   } else if (atomType[0] == 'C') {
2533     if (//Sp3 Carbon
2534       atomType[1] == 'T' ||
2535       strcmp(atomType, "CP1" )==0 ||
2536       strcmp(atomType, "CP2" )==0 ||
2537       strcmp(atomType, "CP3" )==0 ||
2538       strcmp(atomType, "CS"  )==0 ) {
2539       if (numBonds == 1)
2540         return 1;
2541       else if (numBonds == 2)
2542         return 2;
2543       else if (numBonds == 3)
2544         return 3;
2545       else if (numBonds == 4)
2546         return 4;
2547       else
2548         return 1;
2549
2550     } else if (//Sp2
2551       strcmp(atomType, "C"   )==0 ||
2552       strcmp(atomType, "CA"  )==0 ||
2553       strcmp(atomType, "CC"  )==0 ||
2554       strcmp(atomType, "CD"  )==0 ||
2555       strcmp(atomType, "CN"  )==0 ||
2556       strcmp(atomType, "CY"  )==0 ||
2557       strcmp(atomType, "C3"  )==0 ||
2558       strcmp(atomType, "CE1" )==0 ||
2559       strcmp(atomType, "CE2" )==0 ||
2560       strcmp(atomType, "CST" )==0 ||
2561       strcmp(atomType, "CAP" )==0 ||
2562       strcmp(atomType, "COA" )==0 ||
2563       strcmp(atomType, "CPT" )==0 ||
2564       strcmp(atomType, "CPH1")==0 ||
2565       strcmp(atomType, "CPH2")==0
2566       ) {
2567       if (numBonds == 2)
2568         return 5;
2569       else if (numBonds == 3)
2570         return 6;
2571       else
2572         return 1;
2573     } else { // other Carbon
2574         return 1;
2575     }
2576
2577   //Nitrogen
2578   } else if (atomType[0] == 'N') {
2579     if (//Sp3 Nitrogen
2580       //strcmp(atomType, "N"   )==0 ||
2581       //strcmp(atomType, "NH1" )==0 ||
2582       //strcmp(atomType, "NH2" )==0 ||
2583       strcmp(atomType, "NH3" )==0 ||
2584       //strcmp(atomType, "NC2" )==0 ||
2585       //strcmp(atomType, "NY"  )==0 ||
2586       strcmp(atomType, "NP"  )==0
2587       ) {
2588       if (numBonds == 1)
2589         return 11;
2590       else if (numBonds == 2)
2591         return 12;
2592       else if (numBonds == 3)
2593         return 13;
2594       else
2595         return 11;
2596
2597     } else if (//SP2 Nitrogen
2598       strcmp(atomType, "NY"  )==0 || //
2599       strcmp(atomType, "NC2" )==0 || //
2600       strcmp(atomType, "N"   )==0 || //
2601       strcmp(atomType, "NH1" )==0 || //
2602       strcmp(atomType, "NH2" )==0 || //
2603       strcmp(atomType, "NR1" )==0 ||
2604       strcmp(atomType, "NR2" )==0 ||
2605       strcmp(atomType, "NR3" )==0 ||
2606       strcmp(atomType, "NPH" )==0 ||
2607       strcmp(atomType, "NC"  )==0
2608       ) {
2609       if (numBonds == 1)
2610         return 14;
2611       else if (numBonds == 2)
2612         return 15;
2613       else if (numBonds == 3)
2614         return 16;
2615       else
2616         return 11;
2617     } else { // other Nitrogen
2618       return 11;
2619     }
2620
2621   //Oxygen
2622   } else if (atomType[0] == 'O') {
2623     if (//Sp3 Oxygen
2624       strcmp(atomType, "OH1" )==0 ||
2625       strcmp(atomType, "OS"  )==0 ||
2626       strcmp(atomType, "OC"  )==0 || //
2627       strcmp(atomType, "OT"  )==0
2628       ) {
2629       if (numBonds == 1)
2630         return 7;
2631       else if (numBonds == 2)
2632         return 8;
2633       else
2634         return 7;
2635     } else if ( // Sp2 Oxygen
2636       strcmp(atomType, "O"   )==0 ||
2637       strcmp(atomType, "OB"  )==0 ||
2638       strcmp(atomType, "OST" )==0 ||
2639       strcmp(atomType, "OCA" )==0 ||
2640       strcmp(atomType, "OM"  )==0
2641       ) {
2642       return 9;
2643     } else if ( // SP1 Oxygen
2644       strcmp(atomType, "OC"  )==0
2645       ) {
2646       return 10;
2647     } else { // other Oxygen
2648       return 7;
2649     }
2650
2651   //Sulfur
2652   } else if (atomType[0] == 'S') {
2653       if (numBonds == 1)
2654         return 17;
2655       else
2656         return 18;
2657
2658   //Phosphorus
2659   } else if (atomType[0] == 'P') {
2660       if (numBonds == 3)
2661         return 19;
2662       else if (numBonds == 4)
2663         return 20;
2664       else
2665         return 19;
2666   } else { // unknown atom type
2667     return 5;
2668   }
2669   return 5;
2670 } // getLCPOTypeCharmm
2671
2672 //input type is Charmm/Amber/other
2673 //0 - Charmm/Xplor
2674 //1 - Amber
2675 //2 - Plugin
2676 //3 - Gromacs
2677 void Molecule::assignLCPOTypes(int inputType) {
2678   int *heavyBonds = new int[numAtoms];
2679   for (int i = 0; i < numAtoms; i++)
2680     heavyBonds[i] = 0;
2681   for (int i = 0; i < numBonds; i++ ) {
2682     Bond curBond = bonds[i];
2683     int a1 = bonds[i].atom1;
2684     int a2 = bonds[i].atom2;
2685     if (atoms[a1].mass > 2.f && atoms[a2].mass > 2.f) {
2686       heavyBonds[a1]++;
2687       heavyBonds[a2]++;
2688     }
2689   }
2690
2691   lcpoParamType = new int[numAtoms];
2692
2693   int warning = 0;
2694   for (int i = 0; i < numAtoms; i++) {
2695     //print vdw_type and numbonds
2696
2697     if (inputType == 1) { // Amber
2698       lcpoParamType[i] = getLCPOTypeAmber(atomNames[i].atomtype, heavyBonds[i]);
2699     } else { // Charmm
2700       lcpoParamType[i] = getLCPOTypeCharmm(atomNames[i].atomtype, heavyBonds[i]);
2701     }
2702 /*
2703     CkPrintf("%d MOL: ATOM[%05d] = { %4s %d } : %d\n",
2704       inputType,
2705       i+1,
2706       atomNames[i].atomtype,
2707       heavyBonds[i],
2708       lcpoParamType[i]
2709       );
2710 */
2711     if ( atoms[i].mass < 1.5 && lcpoParamType[i] != 0 ) {
2712       if (atoms[i].status & LonepairAtom) {
2713         warning |= LonepairAtom;
2714         lcpoParamType[i] = 0;  // reset LCPO type for LP
2715       }
2716       else if (atoms[i].status & DrudeAtom) {
2717         warning |= DrudeAtom;
2718         lcpoParamType[i] = 0;  // reset LCPO type for Drude
2719       }
2720       else {
2721         CkPrintf("ERROR in Molecule::assignLCPOTypes(): "
2722             "Light atom given heavy atom LCPO type.\n");
2723       }
2724     }
2725
2726     //CkPrintf("VDW_TYPE %02d %4s\n", atoms[i].vdw_type, atomNames[i].atomtype);
2727   } // for atoms
2728
2729   if (warning & LonepairAtom) {
2730     iout << iWARN << "LONE PAIRS TO BE IGNORED BY SASA\n" << endi;
2731   }
2732   if (warning & DrudeAtom) {
2733     iout << iWARN << "DRUDE PARTICLES TO BE IGNORED BY SASA\n" << endi;
2734   }
2735
2736   delete [] heavyBonds;
2737
2738 } // buildLCPOTable
2739
2740 void Molecule::plgLoadAtomBasics(molfile_atom_t *atomarray){
2741     atoms = new Atom[numAtoms];
2742     atomNames = new AtomNameInfo[numAtoms];
2743     if(simParams->genCompressedPsf) {
2744         atomSegResids = new AtomSegResInfo[numAtoms];
2745     }    
2746     hydrogenGroup.resize(0);
2747
2748     ResidueLookupElem *tmpResLookup = resLookup;
2749
2750     for(int i=0; i<numAtoms; i++) {
2751         int reslength = strlen(atomarray[i].resname)+1;
2752         int namelength = strlen(atomarray[i].name)+1;
2753         int typelength = strlen(atomarray[i].type)+1;
2754         atomNames[i].resname = nameArena->getNewArray(reslength);
2755         atomNames[i].atomname = nameArena->getNewArray(namelength);
2756         atomNames[i].atomtype = nameArena->getNewArray(typelength);
2757         strcpy(atomNames[i].resname, atomarray[i].resname);
2758         strcpy(atomNames[i].atomname, atomarray[i].name);
2759         strcpy(atomNames[i].atomtype, atomarray[i].type);
2760
2761         atoms[i].mass = atomarray[i].mass;
2762         atoms[i].charge = atomarray[i].charge;
2763         atoms[i].status = UnknownAtom;
2764
2765         //add this atom to residue lookup table
2766         if(tmpResLookup) {
2767             tmpResLookup = tmpResLookup->append(atomarray[i].segid, atomarray[i].resid, i);
2768         }
2769
2770         if(atomSegResids) { //for compressing molecule information
2771             AtomSegResInfo *one = atomSegResids + i;
2772             memcpy(one->segname, atomarray[i].segid, strlen(atomarray[i].segid)+1);
2773             one->resid = atomarray[i].resid;
2774         }
2775         //Determine the type of the atom
2776         if ( simParams->ignoreMass ) {
2777         }else if(atoms[i].mass <= 0.05) {
2778             atoms[i].status |= LonepairAtom;
2779         }else if(atoms[i].mass < 1.0) {
2780             atoms[i].status |= DrudeAtom;
2781         }else if(atoms[i].mass <= 3.5) {
2782             atoms[i].status |= HydrogenAtom;
2783         }else if((atomNames[i].atomname[0] == 'O') &&
2784                  (atoms[i].mass>=14.0) && (atoms[i].mass<=18.0)){
2785             atoms[i].status |= OxygenAtom;
2786         }
2787         //Look up the vdw constants for this atom
2788         params->assign_vdw_index(atomNames[i].atomtype, &atoms[i]);
2789     }
2790 }
2791
2792 void Molecule::plgLoadBonds(int *from, int *to){
2793     bonds = new Bond[numBonds];
2794     int realNumBonds = 0;
2795     for(int i=0; i<numBonds; i++) {
2796         Bond *thisBond = bonds+realNumBonds;
2797         thisBond->atom1 = from[i]-1;
2798         thisBond->atom2 = to[i]-1;
2799
2800         params->assign_bond_index(
2801             atomNames[thisBond->atom1].atomtype,
2802             atomNames[thisBond->atom2].atomtype,
2803             thisBond);
2804
2805         //Make sure this isn't a fake bond meant for shake in x-plor
2806         Real k, x0;
2807         params->get_bond_params(&k, &x0, thisBond->bond_type);
2808         if (is_lonepairs_psf) {
2809             //need to retain Lonepair bonds for Drude
2810             if(k!=0. || is_lp(thisBond->atom1) || 
2811                is_lp(thisBond->atom2)) {               
2812                 realNumBonds++;
2813             }
2814         }else{
2815             if(k != 0.) realNumBonds++;
2816         }
2817     }
2818
2819     if(numBonds != realNumBonds) {
2820         iout << iWARN << "Ignored" << numBonds-realNumBonds <<
2821             "bonds with zero force constants.\n" <<endi;
2822         iout << iWARN << "Will get H-H distance in rigid H20 from H-O-H angle.\n" <<endi;
2823     }
2824     numBonds = realNumBonds;
2825 }
2826
2827 void Molecule::plgLoadAngles(int *plgAngles)
2828 {    
2829     angles=new Angle[numAngles];
2830     int *atomid = plgAngles;
2831     int numRealAngles = 0;
2832     for(int i=0; i<numAngles; i++) {
2833         Angle *thisAngle = angles+numRealAngles;
2834         thisAngle->atom1 = atomid[0]-1;
2835         thisAngle->atom2 = atomid[1]-1;
2836         thisAngle->atom3 = atomid[2]-1;
2837         atomid += 3;
2838
2839         params->assign_angle_index(
2840             atomNames[thisAngle->atom1].atomtype,
2841             atomNames[thisAngle->atom2].atomtype,
2842             atomNames[thisAngle->atom3].atomtype,
2843                                 thisAngle, simParams->alchOn ? -1 : 0);
2844         if ( thisAngle->angle_type == -1 ) {
2845           iout << iWARN << "ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
2846                << endi;
2847         }
2848
2849         Real k, t0, k_ub, r_ub;
2850         if ( thisAngle->angle_type == -1 ) { k = -1.;  k_ub = -1.; } else
2851         params->get_angle_params(&k, &t0, &k_ub, &r_ub, thisAngle->angle_type);
2852         if(k!=0. || k_ub!=0.) numRealAngles++;
2853     }
2854
2855     if(numAngles != numRealAngles) {
2856         iout << iWARN << "Ignored" << numAngles-numRealAngles << 
2857             " angles with zero force constants.\n" << endi; 
2858     }
2859     numAngles = numRealAngles;
2860 }
2861
2862 void Molecule::plgLoadDihedrals(int *plgDihedrals)
2863 {
2864     std::map< std::string, int > cache;
2865
2866     int lastAtomIds[4];
2867     int multiplicity = 1; //multiplicity of the current bond
2868
2869     lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
2870     dihedrals = new Dihedral[numDihedrals];
2871     int numRealDihedrals = 0;
2872     int *atomid = plgDihedrals;
2873     for(int i=0; i<numDihedrals; i++, atomid+=4) {
2874         Dihedral *thisDihedral = dihedrals + numRealDihedrals;
2875         Bool duplicate_bond = TRUE;
2876         for(int j=0; j<4; j++) {
2877             if(atomid[j] != lastAtomIds[j]) {
2878                 duplicate_bond = FALSE;
2879             }
2880             lastAtomIds[j] = atomid[j];
2881         }
2882
2883         if(duplicate_bond) {
2884             multiplicity++;
2885             if(multiplicity==2) {
2886                 numMultipleDihedrals++;
2887             }
2888         }else{
2889             multiplicity=1;
2890             numRealDihedrals++;
2891         }
2892
2893         thisDihedral->atom1 = atomid[0]-1;
2894         thisDihedral->atom2 = atomid[1]-1;
2895         thisDihedral->atom3 = atomid[2]-1;
2896         thisDihedral->atom4 = atomid[3]-1;
2897
2898       char query[128];
2899       sprintf(query,"%10s :: %10s :: %10s :: %10s :: %d",
2900         atomNames[atomid[0]-1].atomtype,
2901         atomNames[atomid[1]-1].atomtype,
2902         atomNames[atomid[2]-1].atomtype,
2903         atomNames[atomid[3]-1].atomtype,
2904         multiplicity);
2905       auto search = cache.find(query);
2906       if ( search != cache.end() ) { 
2907         thisDihedral->dihedral_type = search->second;
2908       } else {
2909         params->assign_dihedral_index(
2910             atomNames[atomid[0]-1].atomtype,
2911             atomNames[atomid[1]-1].atomtype,
2912             atomNames[atomid[2]-1].atomtype,
2913             atomNames[atomid[3]-1].atomtype,
2914             thisDihedral, multiplicity, simParams->alchOn ? -1 : 0);
2915         if ( thisDihedral->dihedral_type == -1 ) {
2916           iout << iWARN << "ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
2917                << endi;
2918         }
2919         cache[query] = thisDihedral->dihedral_type;
2920       }
2921     }
2922
2923     numDihedrals = numRealDihedrals;
2924 }
2925
2926 void Molecule::plgLoadImpropers(int *plgImpropers)
2927 {
2928     int lastAtomIds[4];
2929     int multiplicity = 1; //multiplicity of the current bond
2930
2931     lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
2932     impropers = new Improper[numImpropers];
2933     int numRealImpropers = 0;
2934     int *atomid = plgImpropers;
2935     for(int i=0; i<numImpropers; i++, atomid+=4) {
2936         Improper *thisImproper = impropers + numRealImpropers;
2937         Bool duplicate_bond = TRUE;
2938         for(int j=0; j<4; j++) {
2939             if(atomid[j] != lastAtomIds[j]) {
2940                 duplicate_bond = FALSE;
2941             }
2942             lastAtomIds[j] = atomid[j];
2943         }
2944
2945         if(duplicate_bond) {
2946             multiplicity++;
2947             if(multiplicity==2) {
2948                 numMultipleImpropers++;
2949             }
2950         }else{
2951             multiplicity=1;
2952             numRealImpropers++;
2953         }
2954
2955         thisImproper->atom1 = atomid[0]-1;
2956         thisImproper->atom2 = atomid[1]-1;
2957         thisImproper->atom3 = atomid[2]-1;
2958         thisImproper->atom4 = atomid[3]-1;
2959
2960         params->assign_improper_index(
2961             atomNames[atomid[0]-1].atomtype,
2962             atomNames[atomid[1]-1].atomtype,
2963             atomNames[atomid[2]-1].atomtype,
2964             atomNames[atomid[3]-1].atomtype,
2965             thisImproper, multiplicity);
2966     }
2967
2968     numImpropers = numRealImpropers;
2969 }
2970
2971 void Molecule::plgLoadCrossterms(int *plgCterms)
2972 {
2973     int lastAtomIds[8];    
2974
2975     for(int i=0; i<8; i++)
2976         lastAtomIds[i]=-1;
2977     
2978     crossterms = new Crossterm[numCrossterms];
2979     int numRealCrossterms = 0;
2980     int *atomid = plgCterms;
2981     for(int i=0; i<numCrossterms; i++, atomid+=8) {
2982         Crossterm *thisCrossterm = crossterms + numRealCrossterms;
2983         Bool duplicate_bond = TRUE;
2984         for(int j=0; j<8; j++) {
2985             if(atomid[j] != lastAtomIds[j]) {
2986                 duplicate_bond = FALSE;
2987             }
2988             lastAtomIds[j] = atomid[j];
2989         }
2990
2991         if(duplicate_bond) {
2992             iout << iWARN <<"Duplicate cross-term detected.\n" << endi;
2993         } else
2994             numRealCrossterms++;
2995
2996         thisCrossterm->atom1 = atomid[0]-1;
2997         thisCrossterm->atom2 = atomid[1]-1;
2998         thisCrossterm->atom3 = atomid[2]-1;
2999         thisCrossterm->atom4 = atomid[3]-1;
3000         thisCrossterm->atom5 = atomid[4]-1;
3001         thisCrossterm->atom6 = atomid[5]-1;
3002         thisCrossterm->atom7 = atomid[6]-1;
3003         thisCrossterm->atom8 = atomid[7]-1;
3004
3005         params->assign_crossterm_index(
3006             atomNames[atomid[0]-1].atomtype,
3007             atomNames[atomid[1]-1].atomtype,
3008             atomNames[atomid[2]-1].atomtype,
3009             atomNames[atomid[3]-1].atomtype,
3010             atomNames[atomid[4]-1].atomtype,
3011             atomNames[atomid[5]-1].atomtype,
3012             atomNames[atomid[6]-1].atomtype,
3013             atomNames[atomid[7]-1].atomtype,
3014             thisCrossterm);
3015     }
3016
3017     numCrossterms = numRealCrossterms;
3018 }
3019
3020 void Molecule::setOccupancyData(molfile_atom_t *atomarray){
3021     occupancy = new float[numAtoms];
3022     for(int i=0; i<numAtoms; i++) {
3023         occupancy[i] = atomarray[i].occupancy;
3024     }
3025 }
3026
3027 void Molecule::setBFactorData(molfile_atom_t *atomarray){
3028     bfactor = new float[numAtoms];
3029     for(int i=0; i<numAtoms; i++) {
3030         bfactor[i] = atomarray[i].bfactor;
3031     }
3032 }
3033
3034     /************************************************************************/
3035     /*                  */
3036     /*      FUNCTION build_lists_by_atom      */
3037     /*                  */
3038     /*  This function builds O(NumAtoms) arrays that store the bonds,   */
3039     /*  angles, dihedrals, and impropers, that each atom is involved in.    */
3040     /*  This is a space hog, but VERY fast.  This will certainly have to    */
3041     /*  change to make things scalable in memory, but for now, speed is the */
3042     /*  thing!                */
3043     /*                  */
3044     /************************************************************************/
3045     void Molecule::build_lists_by_atom()       
3046     {
3047        register int i;      //  Loop counter
3048
3049        register int numFixedAtoms = this->numFixedAtoms;
3050        // if we want forces on fixed atoms then just pretend
3051        // there are none for the purposes of this routine;
3052        if ( simParams->fixedAtomsForces ) numFixedAtoms = 0;
3053
3054 //fepb
3055 //     int numFepInitial = this->numFepInitial;
3056 //     int numFepFinal = this->numFepFinal;
3057 //fepe
3058        tmpArena = new ObjectArena<int32>;
3059        bondsWithAtom = new int32 *[numAtoms];
3060        cluster = new int32 [numAtoms];
3061        clusterSize = new int32 [numAtoms];
3062
3063        bondsByAtom = new int32 *[numAtoms];
3064        anglesByAtom = new int32 *[numAtoms];
3065        dihedralsByAtom = new int32 *[numAtoms];
3066        impropersByAtom = new int32 *[numAtoms];
3067        crosstermsByAtom = new int32 *[numAtoms];
3068
3069        exclusionsByAtom = new int32 *[numAtoms];
3070        fullExclusionsByAtom = new int32 *[numAtoms];
3071        modExclusionsByAtom = new int32 *[numAtoms];
3072
3073        // JLai
3074        gromacsPairByAtom = new int32 *[numAtoms];
3075        // End of JLai
3076
3077        int32 *byAtomSize = new int32[numAtoms];
3078
3079        const int pair_self = 
3080          simParams->pairInteractionOn ? simParams->pairInteractionSelf : 0;
3081
3082        DebugM(3,"Building bond lists.\n");
3083     
3084        //  Build the bond lists
3085        for (i=0; i<numAtoms; i++)
3086        {
3087          byAtomSize[i] = 0;
3088        }
3089        for (i=0; i<numRealBonds; i++)
3090        {
3091          byAtomSize[bonds[i].atom1]++;
3092          byAtomSize[bonds[i].atom2]++;
3093        }
3094        for (i=0; i<numAtoms; i++)
3095        {
3096          bondsWithAtom[i] = tmpArena->getNewArray(byAtomSize[i]+1);
3097          bondsWithAtom[i][byAtomSize[i]] = -1;
3098          byAtomSize[i] = 0;
3099        }
3100        for (i=0; i<numRealBonds; i++)
3101        {
3102          int a1 = bonds[i].atom1;
3103          int a2 = bonds[i].atom2;
3104          bondsWithAtom[a1][byAtomSize[a1]++] = i;
3105          bondsWithAtom[a2][byAtomSize[a2]++] = i;
3106        }
3107
3108         
3109        // Updates all bond, angle, dihedral, improper and crossterm
3110        // to reflect the QM region (which can't have any of there terms)
3111        if (simParams->qmForcesOn) {
3112            
3113            DebugM(3,"Calculating exclusions for QM simulation.\n");
3114            build_exclusions();
3115            
3116            delete_qm_bonded() ;
3117            
3118            DebugM(3,"Re-Building bond lists.\n");
3119            
3120            // We re-calculate the bondsWithAtom list for cluster 
3121            // info calculation below.
3122            for (i=0; i<numAtoms; i++)
3123            {
3124              byAtomSize[i] = 0;
3125            }
3126            for (i=0; i<numRealBonds; i++)
3127            {
3128              byAtomSize[bonds[i].atom1]++;
3129              byAtomSize[bonds[i].atom2]++;
3130            }
3131            for (i=0; i<numAtoms; i++)
3132            {
3133              bondsWithAtom[i][byAtomSize[i]] = -1;
3134              byAtomSize[i] = 0;
3135            }
3136            for (i=0; i<numRealBonds; i++)
3137            {
3138              int a1 = bonds[i].atom1;
3139              int a2 = bonds[i].atom2;
3140              bondsWithAtom[a1][byAtomSize[a1]++] = i;
3141              bondsWithAtom[a2][byAtomSize[a2]++] = i;
3142            }
3143        }
3144         
3145        //  Build cluster information (contiguous clusters)
3146        for (i=0; i<numAtoms; i++) {
3147          cluster[i] = i;
3148        }
3149        for (i=0; i<numAtoms; i++) {
3150          int ci = i;
3151          while ( cluster[ci] != ci ) ci = cluster[ci];
3152          for ( int32 *b = bondsWithAtom[i]; *b != -1; ++b ) {
3153            int a = bonds[*b].atom1;
3154            if ( a == i ) a = bonds[*b].atom2;
3155            if ( a > i ) {
3156              int ca = a;
3157              while ( cluster[ca] != ca ) ca = cluster[ca];
3158              if ( ca > ci ) cluster[ca] = cluster[ci];
3159              else cluster[ci] = cluster[ca];
3160            }
3161          }
3162        }
3163        while ( 1 ) {
3164          int allok = 1;
3165          for (i=0; i<numAtoms; i++) {
3166            int ci = cluster[i];
3167            if ( cluster[ci] != ci ) {
3168              allok = 0;
3169              cluster[i] = cluster[ci];
3170            }
3171          }
3172          if ( allok ) break;
3173        }
3174        
3175        for (i=0; i<numAtoms; i++) {
3176          clusterSize[i] = 0;
3177        }       
3178        for (i=0; i<numAtoms; i++) {           
3179          clusterSize[cluster[i]] += 1;
3180        }
3181
3182 /*
3183        //Getting number of clusters for debugging
3184        int numClusters=0;
3185        for(int i=0; i<numAtoms; i++){
3186            if(clusterSize[i]!=0) numClusters++;
3187        }
3188        printf("Num of clusters: %d\n", numClusters);
3189 */
3190
3191        //  Build the bond lists
3192        for (i=0; i<numAtoms; i++)
3193        {
3194          byAtomSize[i] = 0;
3195        }
3196        numCalcBonds = 0;
3197        for (i=0; i<numBonds; i++)
3198        {
3199          if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3200                             && fixedAtomFlags[bonds[i].atom2] ) continue;
3201    
3202          if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
3203          byAtomSize[bonds[i].atom1]++;
3204          numCalcBonds++;
3205        }
3206        for (i=0; i<numAtoms; i++)
3207        {
3208          bondsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3209          bondsByAtom[i][byAtomSize[i]] = -1;
3210          byAtomSize[i] = 0;
3211        }
3212        for (i=0; i<numBonds; i++)
3213        {
3214          if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3215                             && fixedAtomFlags[bonds[i].atom2] ) continue;
3216          if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
3217          int a1 = bonds[i].atom1;
3218          bondsByAtom[a1][byAtomSize[a1]++] = i;
3219        }
3220        for (i=0; i<numBonds; ++i) {
3221          int a1 = bonds[i].atom1;
3222          int a2 = bonds[i].atom2;
3223          int j;
3224          if ( a1 == a2 ) {
3225            char buff[512];
3226            sprintf(buff,"Atom %d is bonded to itself", a1+1);
3227            NAMD_die(buff);
3228          }
3229          for ( j = 0; j < byAtomSize[a1]; ++j ) {
3230            int b = bondsByAtom[a1][j];
3231            int ba1 = bonds[b].atom1;
3232            int ba2 = bonds[b].atom2;
3233            if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3234              char buff[512];
3235              sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3236              NAMD_die(buff);
3237            }
3238          }
3239          for ( j = 0; j < byAtomSize[a2]; ++j ) {
3240            int b = bondsByAtom[a2][j];
3241            int ba1 = bonds[b].atom1;
3242            int ba2 = bonds[b].atom2;
3243            if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3244              char buff[512];
3245              sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3246              NAMD_die(buff);
3247            }
3248          }
3249        }
3250        
3251        DebugM(3,"Building angle lists.\n");
3252     
3253        //  Build the angle lists
3254        for (i=0; i<numAtoms; i++)
3255        {
3256          byAtomSize[i] = 0;
3257        }
3258        numCalcAngles = 0;
3259        for (i=0; i<numAngles; i++)
3260        {
3261          if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3262                             && fixedAtomFlags[angles[i].atom2]
3263                             && fixedAtomFlags[angles[i].atom3] ) continue;
3264          if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
3265          byAtomSize[angles[i].atom1]++;
3266          numCalcAngles++;
3267        }
3268        for (i=0; i<numAtoms; i++)
3269        {
3270          anglesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3271          anglesByAtom[i][byAtomSize[i]] = -1;
3272          byAtomSize[i] = 0;
3273        }
3274        for (i=0; i<numAngles; i++)
3275        {
3276          if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
3277          if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3278                             && fixedAtomFlags[angles[i].atom2]
3279                             && fixedAtomFlags[angles[i].atom3] ) continue;
3280          int a1 = angles[i].atom1;
3281          anglesByAtom[a1][byAtomSize[a1]++] = i;
3282        }
3283        
3284        DebugM(3,"Building improper lists.\n");
3285     
3286        //  Build the improper lists
3287        for (i=0; i<numAtoms; i++)
3288        {
3289          byAtomSize[i] = 0;
3290        }
3291        numCalcImpropers = 0;
3292        for (i=0; i<numImpropers; i++)
3293        {
3294          if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3295                             && fixedAtomFlags[impropers[i].atom2]
3296                             && fixedAtomFlags[impropers[i].atom3]
3297                             && fixedAtomFlags[impropers[i].atom4] ) continue;
3298          if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
3299          byAtomSize[impropers[i].atom1]++;
3300          numCalcImpropers++;
3301        }
3302        for (i=0; i<numAtoms; i++)
3303        {
3304          impropersByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3305          impropersByAtom[i][byAtomSize[i]] = -1;
3306          byAtomSize[i] = 0;
3307        }
3308        for (i=0; i<numImpropers; i++)
3309        {
3310          if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3311                             && fixedAtomFlags[impropers[i].atom2]
3312                             && fixedAtomFlags[impropers[i].atom3]
3313                             && fixedAtomFlags[impropers[i].atom4] ) continue;
3314          if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
3315          int a1 = impropers[i].atom1;
3316          impropersByAtom[a1][byAtomSize[a1]++] = i;
3317        }
3318        
3319        DebugM(3,"Building dihedral lists.\n");
3320     
3321        //  Build the dihedral lists
3322        for (i=0; i<numAtoms; i++)
3323        {
3324          byAtomSize[i] = 0;
3325        }
3326        numCalcDihedrals = 0;
3327        for (i=0; i<numDihedrals; i++)
3328        {
3329          if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3330                             && fixedAtomFlags[dihedrals[i].atom2]
3331                             && fixedAtomFlags[dihedrals[i].atom3]
3332                             && fixedAtomFlags[dihedrals[i].atom4] ) continue;
3333          if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
3334          byAtomSize[dihedrals[i].atom1]++;
3335          numCalcDihedrals++;
3336        }
3337        for (i=0; i<numAtoms; i++)
3338        {
3339          dihedralsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3340          dihedralsByAtom[i][byAtomSize[i]] = -1;
3341          byAtomSize[i] = 0;
3342        }
3343        for (i=0; i<numDihedrals; i++)
3344        {
3345          if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3346                             && fixedAtomFlags[dihedrals[i].atom2]
3347                             && fixedAtomFlags[dihedrals[i].atom3]
3348                             && fixedAtomFlags[dihedrals[i].atom4] ) continue;
3349          if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
3350          int a1 = dihedrals[i].atom1;
3351          dihedralsByAtom[a1][byAtomSize[a1]++] = i;
3352        }
3353     
3354        DebugM(3,"Building crossterm lists.\n");
3355     
3356        //  Build the crossterm lists
3357        for (i=0; i<numAtoms; i++)
3358        {
3359          byAtomSize[i] = 0;
3360        }
3361        numCalcCrossterms = 0;
3362        for (i=0; i<numCrossterms; i++)
3363        {
3364          if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3365                             && fixedAtomFlags[crossterms[i].atom2]
3366                             && fixedAtomFlags[crossterms[i].atom3]
3367                             && fixedAtomFlags[crossterms[i].atom4]
3368                             && fixedAtomFlags[crossterms[i].atom5]
3369                             && fixedAtomFlags[crossterms[i].atom6]
3370                             && fixedAtomFlags[crossterms[i].atom7]
3371                             && fixedAtomFlags[crossterms[i].atom8] ) continue;
3372          if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
3373          byAtomSize[crossterms[i].atom1]++;
3374          numCalcCrossterms++;
3375        }
3376        for (i=0; i<numAtoms; i++)
3377        {
3378          crosstermsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3379          crosstermsByAtom[i][byAtomSize[i]] = -1;
3380          byAtomSize[i] = 0;
3381        }
3382        for (i=0; i<numCrossterms; i++)
3383        {
3384          if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3385                             && fixedAtomFlags[crossterms[i].atom2]
3386                             && fixedAtomFlags[crossterms[i].atom3]
3387                             && fixedAtomFlags[crossterms[i].atom4]
3388                             && fixedAtomFlags[crossterms[i].atom5]
3389                             && fixedAtomFlags[crossterms[i].atom6]
3390                             && fixedAtomFlags[crossterms[i].atom7]
3391                             && fixedAtomFlags[crossterms[i].atom8] ) continue;
3392          if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
3393          int a1 = crossterms[i].atom1;
3394          crosstermsByAtom[a1][byAtomSize[a1]++] = i;
3395        }
3396
3397        // DRUDE: init lphostIndexes array
3398        if (is_lonepairs_psf) {
3399          // allocate lone pair host index array only if we need it!
3400          DebugM(3,"Initializing lone pair host index array.\n");
3401          lphostIndexes = new int32[numAtoms];
3402          for (i = 0;  i < numAtoms;  i++) {
3403            lphostIndexes[i] = -1;
3404          }
3405          for (i = 0;  i < numLphosts;  i++) {
3406            int32 index = lphosts[i].atom1;
3407            lphostIndexes[index] = i;
3408          }
3409        }
3410        // DRUDE
3411     
3412        // JLai
3413        DebugM(3,"Building gromacsPair lists.\n");
3414     
3415        //  Build the gromacsPair lists
3416        for (i=0; i<numAtoms; i++)
3417        {
3418          byAtomSize[i] = 0;
3419        }
3420        numCalcLJPair = 0;
3421        for (i=0; i<numLJPair; i++)
3422        {
3423          if ( numFixedAtoms && fixedAtomFlags[gromacsPair[i].atom1]
3424                             && fixedAtomFlags[gromacsPair[i].atom2] ) continue;
3425          if ( pair_self && fepAtomFlags[gromacsPair[i].atom1] != 1) continue;
3426          byAtomSize[gromacsPair[i].atom1]++;
3427          numCalcLJPair++;
3428        }
3429        for (i=0; i<numAtoms; i++)
3430        {
3431          gromacsPairByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3432          gromacsPairByAtom[i][byAtomSize[i]] = -1;
3433          byAtomSize[i] = 0;
3434        }
3435        for (i=0; i<numLJPair; i++)
3436        {
3437            if ( numFixedAtoms && fixedAtomFlags[gromacsPair[i].atom1]
3438                 && fixedAtomFlags[gromacsPair[i].atom2] ) continue;
3439            if ( pair_self && fepAtomFlags[gromacsPair[i].atom1] != 1) continue;
3440            int a1 = gromacsPair[i].atom1;
3441            gromacsPairByAtom[a1][byAtomSize[a1]++] = i;
3442            }
3443
3444        // End of JLai
3445
3446        DebugM(3,"Building exclusion data.\n");
3447     
3448        //  Build the arrays of exclusions for each atom
3449        if (! simParams->qmForcesOn)
3450        build_exclusions();
3451
3452        //  Remove temporary structures
3453        delete [] bondsWithAtom;  bondsWithAtom = 0;
3454        delete tmpArena;  tmpArena = 0;
3455
3456        if (exclusions != NULL)
3457       delete [] exclusions;
3458
3459        // 1-4 exclusions which are also fully excluded were eliminated by hash table
3460        numTotalExclusions = exclusionSet.size();
3461        if ( ! CkMyPe() ) {
3462          iout << iINFO << "ADDED " << (numTotalExclusions - numExclusions) << " IMPLICIT EXCLUSIONS\n" << endi;
3463        }
3464        exclusions = new Exclusion[numTotalExclusions];
3465        UniqueSetIter<Exclusion> exclIter(exclusionSet);
3466        for ( exclIter=exclIter.begin(),i=0; exclIter != exclIter.end(); exclIter++,i++ )
3467        {
3468          exclusions[i] = *exclIter;
3469        }
3470        // Free exclusionSet storage
3471        // exclusionSet.clear(1);
3472        exclusionSet.clear();
3473
3474        DebugM(3,"Building exclusion lists.\n");
3475     
3476        for (i=0; i<numAtoms; i++)
3477        {
3478          byAtomSize[i] = 0;
3479        }
3480        numCalcExclusions = 0;
3481        numCalcFullExclusions = 0;
3482        for (i=0; i<numTotalExclusions; i++)
3483        {
3484          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3485                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
3486          byAtomSize[exclusions[i].atom1]++;
3487          numCalcExclusions++;
3488          if ( ! exclusions[i].modified ) numCalcFullExclusions++;
3489        }
3490
3491        for (i=0; i<numAtoms; i++)
3492        {
3493          exclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3494          exclusionsByAtom[i][byAtomSize[i]] = -1;
3495          byAtomSize[i] = 0;
3496        }
3497        for (i=0; i<numTotalExclusions; i++)
3498        {
3499          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3500                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
3501          int a1 = exclusions[i].atom1;
3502          exclusionsByAtom[a1][byAtomSize[a1]++] = i;
3503        }
3504
3505        int32 *byAtomSize2 = new int32[numAtoms];
3506
3507        for (i=0; i<numAtoms; i++)
3508        {
3509          byAtomSize[i] = 0;
3510          byAtomSize2[i] = 0;
3511        }
3512
3513        for (i=0; i<numTotalExclusions; i++)
3514        {
3515          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3516                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
3517          if ( exclusions[i].modified ) {
3518            byAtomSize2[exclusions[i].atom1]++;
3519            byAtomSize2[exclusions[i].atom2]++;
3520          } else {
3521            byAtomSize[exclusions[i].atom1]++;
3522            byAtomSize[exclusions[i].atom2]++;
3523          }
3524        }
3525
3526        for (i=0; i<numAtoms; i++)
3527        {
3528          fullExclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3529          fullExclusionsByAtom[i][0] = 0;
3530          modExclusionsByAtom[i] = arena->getNewArray(byAtomSize2[i]+1);
3531          modExclusionsByAtom[i][0] = 0;
3532        }
3533
3534        for (i=0; i<numTotalExclusions; i++)
3535        {
3536          int a1 = exclusions[i].atom1;
3537          int a2 = exclusions[i].atom2;
3538          if ( numFixedAtoms && fixedAtomFlags[a1]
3539                             && fixedAtomFlags[a2] ) continue;
3540          int32 *l1, *l2;
3541          if ( exclusions[i].modified ) {
3542            l1 = modExclusionsByAtom[a1];
3543            l2 = modExclusionsByAtom[a2];
3544          } else {
3545            l1 = fullExclusionsByAtom[a1];
3546            l2 = fullExclusionsByAtom[a2];
3547          }
3548          l1[++(*l1)] = a2;
3549          l2[++(*l2)] = a1;
3550        }
3551
3552        if ( ! CkMyPe() && simParams->printExclusions ) {
3553          for (i=0; i<numAtoms; i++) {
3554            int32 *lf = fullExclusionsByAtom[i];
3555            iout << "EXCL " << i << " FULL";
3556            int nf = *(lf++);
3557            for ( int j = 0; j < nf; ++j ) {
3558              iout << " " << *(lf++);
3559            }
3560            iout << "\n";
3561            int32 *lm = modExclusionsByAtom[i];
3562            iout << "EXCL " << i << " MOD";
3563            int nm = *(lm++);
3564            for ( int j = 0; j < nm; ++j ) {
3565              iout << " " << *(lm++);
3566            }
3567            iout << "\n" << endi;
3568          }
3569        }
3570
3571        // DRUDE
3572        if (is_drude_psf || simParams->drudeOn) {
3573
3574          // build Thole (screened Coulomb) correction terms;
3575          // they are constructed implicitly from exclusions
3576
3577          // free the previous Thole array if already allocated
3578          if (tholes != NULL) delete[] tholes;
3579          numTholes = 0;
3580
3581          // count the number of Thole terms
3582          for (i = 0;  i < numTotalExclusions;  i++) {
3583            /* skip over the modified exclusions */
3584            if (exclusions[i].modified) continue;
3585            int a1 = exclusions[i].atom1;
3586            int a2 = exclusions[i].atom2;
3587            if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3588              numTholes++;
3589            }
3590          }
3591
3592          // allocate space for Thole terms
3593          if (numTholes != 0) tholes = new Thole[numTholes];
3594          else tholes = NULL;
3595          int nt = 0;
3596
3597          Real c = COULOMB*simParams->nonbondedScaling/simParams->dielectric;
3598
3599          // store Thole terms
3600          for (i = 0;  i < numTotalExclusions;  i++) {
3601            /* skip over the modified exclusions */
3602            if (exclusions[i].modified) continue;
3603            int a1 = exclusions[i].atom1;
3604            int a2 = exclusions[i].atom2;
3605            // exclusions are stored with a1 < a2
3606            if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3607              Real thsum = drudeConsts[a1].thole + drudeConsts[a2].thole;
3608              Real aprod = drudeConsts[a1].alpha * drudeConsts[a2].alpha;
3609              // guard against having alpha==0
3610              Real apower = (aprod <= 0 ? 0 : powf(aprod, -1.f/6));
3611              tholes[nt].atom1 = a1;
3612              tholes[nt].atom2 = a1+1;
3613              tholes[nt].atom3 = a2;
3614              tholes[nt].atom4 = a2+1;
3615              tholes[nt].aa = apower * thsum;
3616              tholes[nt].qq = c * atoms[a1+1].charge * atoms[a2+1].charge;
3617              nt++;
3618            }
3619          }
3620
3621          // build Thole lists by atom
3622          DebugM(3, "Building Thole correction term lists.\n");
3623          tholesByAtom = new int32 *[numAtoms];
3624
3625          for (i = 0;  i < numAtoms;  i++) {
3626            byAtomSize[i] = 0;
3627          }
3628          numCalcTholes = 0;
3629          for (i = 0;  i < numTholes;  i++) {
3630            if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3631                               && fixedAtomFlags[tholes[i].atom2]
3632                               && fixedAtomFlags[tholes[i].atom3]
3633                               && fixedAtomFlags[tholes[i].atom4] ) continue;
3634            if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1) continue;
3635            byAtomSize[tholes[i].atom1]++;
3636            numCalcTholes++;
3637          }
3638          for (i = 0;  i < numAtoms;  i++) {
3639            tholesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3640            tholesByAtom[i][byAtomSize[i]] = -1;
3641            byAtomSize[i] = 0;
3642          }
3643          for (i = 0;  i < numTholes;  i++) {
3644            if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3645                               && fixedAtomFlags[tholes[i].atom2]
3646                               && fixedAtomFlags[tholes[i].atom3]
3647                               && fixedAtomFlags[tholes[i].atom4] ) continue;
3648            if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1) continue;
3649            int a1 = tholes[i].atom1;
3650            tholesByAtom[a1][byAtomSize[a1]++] = i;
3651          }
3652
3653          // build anisotropic lists by atom
3654          DebugM(3, "Building anisotropic term lists.\n");
3655          anisosByAtom = new int32 *[numAtoms];
3656
3657          for (i = 0;  i < numAtoms;  i++) {
3658            byAtomSize[i] = 0;
3659          }
3660          numCalcAnisos = 0;
3661          for (i = 0;  i < numAnisos;  i++) {
3662            if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
3663                               && fixedAtomFlags[anisos[i].atom2]
3664                               && fixedAtomFlags[anisos[i].atom3]
3665                               && fixedAtomFlags[anisos[i].atom4] ) continue;
3666            if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1) continue;
3667            byAtomSize[anisos[i].atom1]++;
3668