01c4867a7c8d7917809c615be2fe68e93a9c1584
[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   char atom1name[11];  // Atom type for atom #1
1484   char atom2name[11];  // Atom type for atom #2
1485   register int j;      // Loop counter
1486   int num_read=0;    // Number of bonds read so far
1487   int origNumBonds = numBonds;   // number of bonds in file header
1488
1489   /*  Allocate the array to hold the bonds      */
1490   bonds=new Bond[numBonds];
1491
1492   if (bonds == NULL)
1493   {
1494     NAMD_die("memory allocations failed in Molecule::read_bonds");
1495   }
1496
1497   /*  Loop through and read in all the bonds      */
1498   while (num_read < numBonds)
1499   {
1500     /*  Loop and read in the two atom indexes    */
1501     for (j=0; j<2; j++)
1502     {
1503       /*  Read the atom number from the file.         */
1504       /*  Subtract 1 to convert the index from the    */
1505       /*  1 to NumAtoms used in the file to the       */
1506       /*  0 to NumAtoms-1 that we need    */
1507       atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
1508
1509       /*  Check to make sure the index isn't too big  */
1510       if (atom_nums[j] >= numAtoms)
1511       {
1512         char err_msg[128];
1513
1514         sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1515         NAMD_die(err_msg);
1516       }
1517     }
1518
1519     /*  Get the atom type for the two atoms.  When we query */
1520     /*  the parameter object, we need to send the atom type */
1521     /*  that is alphabetically first as atom 1.    */
1522     if (strcasecmp(atomNames[atom_nums[0]].atomtype, 
1523          atomNames[atom_nums[1]].atomtype) < 0)
1524     {
1525       strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
1526       strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
1527     }
1528     else
1529     {
1530       strcpy(atom2name, atomNames[atom_nums[0]].atomtype);
1531       strcpy(atom1name, atomNames[atom_nums[1]].atomtype);
1532     }
1533
1534     /*  Assign the atom indexes to the array element  */
1535     Bond *b = &(bonds[num_read]);
1536     b->atom1=atom_nums[0];
1537     b->atom2=atom_nums[1];
1538
1539     /*  Query the parameter object for the constants for    */
1540     /*  this bond            */
1541     params->assign_bond_index(atom1name, atom2name, b);
1542
1543     /*  Make sure this isn't a fake bond meant for shake in x-plor.  */
1544     Real k, x0;
1545     params->get_bond_params(&k,&x0,b->bond_type);
1546     if (is_lonepairs_psf) {
1547       // need to retain Lonepair bonds for Drude
1548       if ( k == 0. && !is_lp(b->atom1) && !is_lp(b->atom2)) --numBonds;
1549       else ++num_read;
1550     }
1551     else {
1552       if ( k == 0. ) --numBonds;  // fake bond
1553       else ++num_read;  // real bond
1554     }
1555   }
1556
1557   /*  Tell user about our subterfuge  */
1558   if ( numBonds != origNumBonds ) {
1559     iout << iWARN << "Ignored " << origNumBonds - numBonds <<
1560             " bonds with zero force constants.\n" << endi;
1561     iout << iWARN <<
1562         "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
1563   }
1564
1565   return;
1566 }
1567 /*      END OF FUNCTION read_bonds      */
1568
1569 /************************************************************************/
1570 /*                  */
1571 /*      FUNCTION read_angles        */
1572 /*                  */
1573 /*   INPUTS:                */
1574 /*  fd - File descriptor for .psf file        */
1575 /*  params - Parameters object to query for parameters    */
1576 /*                  */
1577 /*  read_angles reads the angle parameters from the .psf file.      */
1578 /*   This section of the .psf file consists of a list of triplets of    */
1579 /*   atom indexes.  Each triplet represents three atoms connected via   */
1580 /*   an angle bond.  The parameter object is queried to obtain the      */
1581 /*   constants for each bond.            */
1582 /*                  */
1583 /************************************************************************/
1584
1585 void Molecule::read_angles(FILE *fd, Parameters *params)
1586
1587 {
1588   int atom_nums[3];  //  Atom numbers for the three atoms
1589   char atom1name[11];  //  Atom type for atom 1
1590   char atom2name[11];  //  Atom type for atom 2
1591   char atom3name[11];  //  Atom type for atom 3
1592   register int j;      //  Loop counter
1593   int num_read=0;    //  Number of angles read so far
1594   int origNumAngles = numAngles;  // Number of angles in file
1595   /*  Alloc the array of angles          */
1596   angles=new Angle[numAngles];
1597
1598   if (angles == NULL)
1599   {
1600     NAMD_die("memory allocation failed in Molecule::read_angles");
1601   }
1602
1603   /*  Loop through and read all the angles      */
1604   while (num_read < numAngles)
1605   {
1606     /*  Loop through the 3 atom indexes in the current angle*/
1607     for (j=0; j<3; j++)
1608     {
1609       /*  Read the atom number from the file.         */
1610       /*  Subtract 1 to convert the index from the    */
1611       /*  1 to NumAtoms used in the file to the       */
1612       /*  0 to NumAtoms-1 that we need    */
1613       atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
1614
1615       /*  Check to make sure the atom index doesn't   */
1616       /*  exceed the Number of Atoms      */
1617       if (atom_nums[j] >= numAtoms)
1618       {
1619         char err_msg[128];
1620
1621         sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1622         NAMD_die(err_msg);
1623       }
1624     }
1625
1626     /*  Place the bond name that is alphabetically first  */
1627     /*  in the atom1name.  This is OK since the order of    */
1628     /*  atom1 and atom3 are interchangable.  And to search  */
1629     /*  the tree of angle parameters, we need the order     */
1630     /*  to be predictable.          */
1631     if (strcasecmp(atomNames[atom_nums[0]].atomtype, 
1632          atomNames[atom_nums[2]].atomtype) < 0)
1633     {
1634       strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
1635       strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
1636       strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
1637     }
1638     else
1639     {
1640       strcpy(atom1name, atomNames[atom_nums[2]].atomtype);
1641       strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
1642       strcpy(atom3name, atomNames[atom_nums[0]].atomtype);
1643     }
1644
1645     /*  Assign the three atom indices      */
1646     angles[num_read].atom1=atom_nums[0];
1647     angles[num_read].atom2=atom_nums[1];
1648     angles[num_read].atom3=atom_nums[2];
1649
1650     /*  Get the constant values for this bond from the  */
1651     /*  parameter object          */
1652     params->assign_angle_index(atom1name, atom2name, 
1653        atom3name, &(angles[num_read]), simParams->alchOn ? -1 : 0);
1654     if ( angles[num_read].angle_type == -1 ) {
1655       iout << iWARN << "ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
1656            << endi;
1657     }
1658
1659     /*  Make sure this isn't a fake angle meant for shake in x-plor.  */
1660     Real k, t0, k_ub, r_ub;
1661     if ( angles[num_read].angle_type == -1 ) { k = -1.;  k_ub = -1.; } else
1662     params->get_angle_params(&k,&t0,&k_ub,&r_ub,angles[num_read].angle_type);
1663     if ( k == 0. && k_ub == 0. ) --numAngles;  // fake angle
1664     else ++num_read;  // real angle
1665   }
1666
1667   /*  Tell user about our subterfuge  */
1668   if ( numAngles != origNumAngles ) {
1669     iout << iWARN << "Ignored " << origNumAngles - numAngles <<
1670             " angles with zero force constants.\n" << endi;
1671   }
1672
1673   return;
1674 }
1675 /*      END OF FUNCTION read_angles      */
1676
1677 /************************************************************************/
1678 /*                  */
1679 /*        FUNCTION read_dihedrals      */
1680 /*                  */
1681 /*   INPUTS:                */
1682 /*  fd - file descriptor for the .psf file        */
1683 /*  params - pointer to parameter object        */
1684 /*                  */
1685 /*  read_dihedreals reads the dihedral section of the .psf file.    */
1686 /*   This section of the file contains a list of quartets of atom       */
1687 /*   numbers.  Each quartet represents a group of atoms that form a     */
1688 /*   dihedral bond.              */
1689 /*                  */
1690 /************************************************************************/
1691
1692 void Molecule::read_dihedrals(FILE *fd, Parameters *params)
1693 {
1694   int atom_nums[4];  // The 4 atom indexes
1695   int last_atom_nums[4];  // Atom numbers from previous bond
1696   char atom1name[11];  // Atom type for atom 1
1697   char atom2name[11];  // Atom type for atom 2
1698   char atom3name[11];  // Atom type for atom 3
1699   char atom4name[11];  // Atom type for atom 4
1700   register int j;      // loop counter
1701   int num_read=0;    // number of dihedrals read so far
1702   int multiplicity=1;  // multiplicity of the current bond
1703   Bool duplicate_bond;  // Is this a duplicate of the last bond
1704   int num_unique=0;   // Number of unique dihedral bonds
1705
1706   //  Initialize the array used to check for duplicate dihedrals
1707   for (j=0; j<4; j++)
1708     last_atom_nums[j] = -1;
1709
1710   /*  Allocate an array to hold the Dihedrals      */
1711   dihedrals = new Dihedral[numDihedrals];
1712
1713   if (dihedrals == NULL)
1714   {
1715     NAMD_die("memory allocation failed in Molecule::read_dihedrals");
1716   }
1717
1718   /*  Loop through and read all the dihedrals      */
1719   while (num_read < numDihedrals)
1720   {
1721     duplicate_bond = TRUE;
1722
1723     /*  Loop through and read the 4 indexes for this bond   */
1724     for (j=0; j<4; j++)
1725     {
1726       /*  Read the atom number from the file.         */
1727       /*  Subtract 1 to convert the index from the    */
1728       /*  1 to NumAtoms used in the file to the       */
1729       /*  0 to NumAtoms-1 that we need    */
1730       atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
1731
1732       /*  Check for an atom index that is too large  */
1733       if (atom_nums[j] >= numAtoms)
1734       {
1735         char err_msg[128];
1736
1737         sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1738         NAMD_die(err_msg);
1739       }
1740
1741       //  Check to see if this atom matches the last bond
1742       if (atom_nums[j] != last_atom_nums[j])
1743       {
1744         duplicate_bond = FALSE;
1745       }
1746
1747       last_atom_nums[j] = atom_nums[j];
1748     }
1749
1750     /*  Get the atom types for the 4 atoms so we can look  */
1751     /*  up the constants in the parameter object    */
1752     strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
1753     strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
1754     strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
1755     strcpy(atom4name, atomNames[atom_nums[3]].atomtype);
1756
1757     //  Check to see if this is really a new bond or just
1758     //  a repeat of the last one
1759     if (duplicate_bond)
1760     {
1761       //  This is a duplicate, so increase the multiplicity
1762       multiplicity++;
1763
1764       if (multiplicity == 2)
1765       {
1766         numMultipleDihedrals++;
1767       }
1768     }
1769     else
1770     {
1771       multiplicity=1;
1772       num_unique++;
1773     }
1774
1775     /*  Assign the atom indexes        */
1776     dihedrals[num_unique-1].atom1=atom_nums[0];
1777     dihedrals[num_unique-1].atom2=atom_nums[1];
1778     dihedrals[num_unique-1].atom3=atom_nums[2];
1779     dihedrals[num_unique-1].atom4=atom_nums[3];
1780
1781     /*  Get the constants for this dihedral bond    */
1782     params->assign_dihedral_index(atom1name, atom2name, 
1783        atom3name, atom4name, &(dihedrals[num_unique-1]),
1784        multiplicity, simParams->alchOn ? -1 : 0);
1785     if ( dihedrals[num_unique-1].dihedral_type == -1 ) {
1786       iout << iWARN << "ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
1787            << endi;
1788     }
1789
1790     num_read++;
1791   }
1792
1793   numDihedrals = num_unique;
1794
1795   return;
1796 }
1797 /*      END OF FUNCTION read_dihedral      */
1798
1799 /************************************************************************/
1800 /*                  */
1801 /*        FUNCTION read_impropers      */
1802 /*                  */
1803 /*   INPUTS:                */
1804 /*  fd - file descriptor for .psf file        */
1805 /*  params - parameter object          */
1806 /*                  */
1807 /*  read_impropers reads the improper section of the .psf file.  */
1808 /*   This section is identical to the dihedral section in that it is    */
1809 /*   made up of a list of quartets of atom indexes that define the      */
1810 /*   atoms that are bonded together.          */
1811 /*                  */
1812 /************************************************************************/
1813
1814 void Molecule::read_impropers(FILE *fd, Parameters *params)
1815 {
1816   int atom_nums[4];  //  Atom indexes for the 4 atoms
1817   int last_atom_nums[4];  //  Atom indexes from previous bond
1818   char atom1name[11];  //  Atom type for atom 1
1819   char atom2name[11];  //  Atom type for atom 2
1820   char atom3name[11];  //  Atom type for atom 3
1821   char atom4name[11];  //  Atom type for atom 4
1822   register int j;      //  Loop counter
1823   int num_read=0;    //  Number of impropers read so far
1824   int multiplicity=1;  // multiplicity of the current bond
1825   Bool duplicate_bond;  // Is this a duplicate of the last bond
1826   int num_unique=0;   // Number of unique dihedral bonds
1827
1828   //  Initialize the array used to look for duplicate improper
1829   //  entries.  Set them all to -1 so we know nothing will match
1830   for (j=0; j<4; j++)
1831     last_atom_nums[j] = -1;
1832
1833   /*  Allocate the array to hold the impropers      */
1834   impropers=new Improper[numImpropers];
1835
1836   if (impropers == NULL)
1837   {
1838     NAMD_die("memory allocation failed in Molecule::read_impropers");
1839   }
1840
1841   /*  Loop through and read all the impropers      */
1842   while (num_read < numImpropers)
1843   {
1844     duplicate_bond = TRUE;
1845
1846     /*  Loop through the 4 indexes for this improper  */
1847     for (j=0; j<4; j++)
1848     {
1849       /*  Read the atom number from the file.         */
1850       /*  Subtract 1 to convert the index from the    */
1851       /*  1 to NumAtoms used in the file to the       */
1852       /*  0 to NumAtoms-1 that we need    */
1853       atom_nums[j]=NAMD_read_int(fd, "IMPROPERS")-1;
1854
1855       /*  Check to make sure the index isn't too big  */
1856       if (atom_nums[j] >= numAtoms)
1857       {
1858         char err_msg[128];
1859
1860         sprintf(err_msg, "IMPROPERS INDEX %d GREATER THAN NATOM %d IN IMPROPERS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
1861         NAMD_die(err_msg);
1862       }
1863
1864       if (atom_nums[j] != last_atom_nums[j])
1865       {
1866         duplicate_bond = FALSE;
1867       }
1868
1869       last_atom_nums[j] = atom_nums[j];
1870     }
1871
1872     /*  Get the atom types so we can look up the parameters */
1873     strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
1874     strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
1875     strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
1876     strcpy(atom4name, atomNames[atom_nums[3]].atomtype);
1877
1878     //  Check to see if this is a duplicate improper
1879     if (duplicate_bond)
1880     {
1881       //  This is a duplicate improper.  So we don't
1882       //  really count this entry, we just update
1883       //  the parameters object
1884       multiplicity++;
1885
1886       if (multiplicity == 2)
1887       {
1888         //  Count the number of multiples.
1889         numMultipleImpropers++;
1890       }
1891     }
1892     else
1893     {
1894       //  Not a duplicate
1895       multiplicity = 1;
1896       num_unique++;
1897     }
1898
1899     /*  Assign the atom indexes        */
1900     impropers[num_unique-1].atom1=atom_nums[0];
1901     impropers[num_unique-1].atom2=atom_nums[1];
1902     impropers[num_unique-1].atom3=atom_nums[2];
1903     impropers[num_unique-1].atom4=atom_nums[3];
1904
1905     /*  Look up the constants for this bond      */
1906     params->assign_improper_index(atom1name, atom2name, 
1907        atom3name, atom4name, &(impropers[num_unique-1]),
1908        multiplicity);
1909
1910     num_read++;
1911   }
1912
1913   //  Now reset the numImpropers value to the number of UNIQUE
1914   //  impropers.  Sure, we waste a few entries in the improper_array
1915   //  on the master node, but it is very little space . . .
1916   numImpropers = num_unique;
1917
1918   return;
1919 }
1920 /*      END OF FUNCTION read_impropers      */
1921
1922 /************************************************************************/
1923 /*                  */
1924 /*        FUNCTION read_crossterms      */
1925 /*                  */
1926 /*   INPUTS:                */
1927 /*  fd - file descriptor for .psf file        */
1928 /*  params - parameter object          */
1929 /*                  */
1930 /*   This section is identical to the dihedral section in that it is    */
1931 /*   made up of a list of quartets of atom indexes that define the      */
1932 /*   atoms that are bonded together.          */
1933 /*                  */
1934 /************************************************************************/
1935
1936 void Molecule::read_crossterms(FILE *fd, Parameters *params)
1937
1938 {
1939   int atom_nums[8];  //  Atom indexes for the 4 atoms
1940   int last_atom_nums[8];  //  Atom indexes from previous bond
1941   char atom1name[11];  //  Atom type for atom 1
1942   char atom2name[11];  //  Atom type for atom 2
1943   char atom3name[11];  //  Atom type for atom 3
1944   char atom4name[11];  //  Atom type for atom 4
1945   char atom5name[11];  //  Atom type for atom 5
1946   char atom6name[11];  //  Atom type for atom 6
1947   char atom7name[11];  //  Atom type for atom 7
1948   char atom8name[11];  //  Atom type for atom 8
1949   register int j;      //  Loop counter
1950   int num_read=0;    //  Number of items read so far
1951   Bool duplicate_bond;  // Is this a duplicate of the last bond
1952
1953   //  Initialize the array used to look for duplicate crossterm
1954   //  entries.  Set them all to -1 so we know nothing will match
1955   for (j=0; j<8; j++)
1956     last_atom_nums[j] = -1;
1957
1958   /*  Allocate the array to hold the cross-terms */
1959   crossterms=new Crossterm[numCrossterms];
1960
1961   if (crossterms == NULL)
1962   {
1963     NAMD_die("memory allocation failed in Molecule::read_crossterms");
1964   }
1965
1966   /*  Loop through and read all the cross-terms      */
1967   while (num_read < numCrossterms)
1968   {
1969     duplicate_bond = TRUE;
1970
1971     /*  Loop through the 8 indexes for this cross-term */
1972     for (j=0; j<8; j++)
1973     {
1974       /*  Read the atom number from the file.         */
1975       /*  Subtract 1 to convert the index from the    */
1976       /*  1 to NumAtoms used in the file to the       */
1977       /*  0 to NumAtoms-1 that we need    */
1978       atom_nums[j]=NAMD_read_int(fd, "CROSS-TERMS")-1;
1979
1980       /*  Check to make sure the index isn't too big  */
1981       if (atom_nums[j] >= numAtoms)
1982       {
1983         char err_msg[128];
1984
1985         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);
1986         NAMD_die(err_msg);
1987       }
1988
1989       if (atom_nums[j] != last_atom_nums[j])
1990       {
1991         duplicate_bond = FALSE;
1992       }
1993
1994       last_atom_nums[j] = atom_nums[j];
1995     }
1996
1997     /*  Get the atom types so we can look up the parameters */
1998     strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
1999     strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
2000     strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
2001     strcpy(atom4name, atomNames[atom_nums[3]].atomtype);
2002     strcpy(atom5name, atomNames[atom_nums[4]].atomtype);
2003     strcpy(atom6name, atomNames[atom_nums[5]].atomtype);
2004     strcpy(atom7name, atomNames[atom_nums[6]].atomtype);
2005     strcpy(atom8name, atomNames[atom_nums[7]].atomtype);
2006
2007     //  Check to see if this is a duplicate term
2008     if (duplicate_bond)
2009     {
2010       iout << iWARN << "Duplicate cross-term detected.\n" << endi;
2011     }
2012
2013     /*  Assign the atom indexes        */
2014     crossterms[num_read].atom1=atom_nums[0];
2015     crossterms[num_read].atom2=atom_nums[1];
2016     crossterms[num_read].atom3=atom_nums[2];
2017     crossterms[num_read].atom4=atom_nums[3];
2018     crossterms[num_read].atom5=atom_nums[4];
2019     crossterms[num_read].atom6=atom_nums[5];
2020     crossterms[num_read].atom7=atom_nums[6];
2021     crossterms[num_read].atom8=atom_nums[7];
2022
2023     /*  Look up the constants for this bond      */
2024     params->assign_crossterm_index(atom1name, atom2name, 
2025        atom3name, atom4name, atom5name, atom6name,
2026        atom7name, atom8name, &(crossterms[num_read]));
2027
2028     if(!duplicate_bond) num_read++;
2029   }
2030
2031   numCrossterms = num_read;
2032
2033   return;
2034 }
2035 /*      END OF FUNCTION read_impropers      */
2036
2037 /************************************************************************/
2038 /*                  */
2039 /*      FUNCTION read_donors        */
2040 /*                  */
2041 /*  read_donors reads in the bond section of the .psf file.  This   */
2042 /*  section contains a list of pairs of numbers where each pair is      */
2043 /*  represents two atoms that are part of an H-bond.  Each atom pair is */
2044 /*  read in.                                                            */
2045 /*                  */
2046 /*  Donor atoms are the heavy atoms to which hydrogens are bonded.      */
2047 /*  There will always be a donor atom for each donor pair.  However,    */
2048 /*  for a united-atom model there may not be an explicit hydrogen       */
2049 /*  present, in which case the second atom index in the pair will be    */
2050 /*  given as 0 in the PSF (and stored as -1 in this program's internal  */
2051 /*  storage).                                                           */
2052 /************************************************************************/
2053
2054 void Molecule::read_donors(FILE *fd)
2055 {
2056   int d[2];               // temporary storage of donor atom index
2057   register int j;      // Loop counter
2058   int num_read=0;    // Number of bonds read so far
2059   int num_no_hydr=0;      // Number of bonds with no hydrogen given
2060
2061   /*  Allocate the array to hold the bonds      */
2062   donors=new Bond[numDonors];
2063
2064   if (donors == NULL)
2065   {
2066     NAMD_die("memory allocations failed in Molecule::read_donors");
2067   }
2068
2069   /*  Loop through and read in all the donors      */
2070   while (num_read < numDonors)
2071   {
2072     /*  Loop and read in the two atom indexes    */
2073     for (j=0; j<2; j++)
2074     {
2075       /*  Read the atom number from the file.         */
2076       /*  Subtract 1 to convert the index from the    */
2077       /*  1 to NumAtoms used in the file to the       */
2078       /*  0 to NumAtoms-1 that we need    */
2079       d[j]=NAMD_read_int(fd, "DONORS")-1;
2080
2081       /*  Check to make sure the index isn't too big  */
2082       if (d[j] >= numAtoms)
2083       {
2084         char err_msg[128];
2085
2086         sprintf(err_msg,
2087     "DONOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE",
2088           d[j]+1, numAtoms, num_read+1);
2089         NAMD_die(err_msg);
2090       }
2091
2092       /*  Check if there is a hydrogen given */
2093       if (d[j] < 0)
2094                           num_no_hydr++;
2095     }
2096
2097     /*  Assign the atom indexes to the array element  */
2098     Bond *b = &(donors[num_read]);
2099     b->atom1=d[0];
2100     b->atom2=d[1];
2101
2102     num_read++;
2103   }
2104
2105   return;
2106 }
2107 /*      END OF FUNCTION read_donors      */
2108
2109
2110 /************************************************************************/
2111 /*                  */
2112 /*      FUNCTION read_acceptors        */
2113 /*                  */
2114 /*  read_acceptors reads in the bond section of the .psf file.      */
2115 /*  This section contains a list of pairs of numbers where each pair is */
2116 /*  represents two atoms that are part of an H-bond.  Each atom pair is */
2117 /*  read in.                                                            */
2118 /*                  */
2119 /*  Acceptor atoms are the heavy atoms to which hydrogens directly      */
2120 /*  orient in a hydrogen bond interaction.  There will always be an     */
2121 /*  acceptor atom for each acceptor pair.  The antecedent atom, to      */
2122 /*  which the acceptor is bound, may not be given in the structure,     */
2123 /*  however, in which case the second atom index in the pair will be    */
2124 /*  given as 0 in the PSF (and stored as -1 in this program's internal  */
2125 /*  storage).                                                           */
2126 /************************************************************************/
2127
2128 void Molecule::read_acceptors(FILE *fd)
2129 {
2130   int d[2];               // temporary storage of atom index
2131   register int j;      // Loop counter
2132   int num_read=0;    // Number of bonds read so far
2133         int num_no_ante=0;      // number of pairs with no antecedent
2134
2135   /*  Allocate the array to hold the bonds      */
2136   acceptors=new Bond[numAcceptors];
2137
2138   if (acceptors == NULL)
2139   {
2140     NAMD_die("memory allocations failed in Molecule::read_acceptors");
2141   }
2142
2143   /*  Loop through and read in all the acceptors      */
2144   while (num_read < numAcceptors)
2145   {
2146     /*  Loop and read in the two atom indexes    */
2147     for (j=0; j<2; j++)
2148     {
2149       /*  Read the atom number from the file.         */
2150       /*  Subtract 1 to convert the index from the    */
2151       /*  1 to NumAtoms used in the file to the       */
2152       /*  0 to NumAtoms-1 that we need    */
2153       d[j]=NAMD_read_int(fd, "ACCEPTORS")-1;
2154
2155       /*  Check to make sure the index isn't too big  */
2156       if (d[j] >= numAtoms)
2157       {
2158         char err_msg[128];
2159
2160         sprintf(err_msg, "ACCEPTOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE", d[j]+1, numAtoms, num_read+1);
2161         NAMD_die(err_msg);
2162       }
2163
2164       /*  Check if there is an antecedent given */
2165       if (d[j] < 0)
2166                           num_no_ante++;
2167     }
2168
2169     /*  Assign the atom indexes to the array element  */
2170     Bond *b = &(acceptors[num_read]);
2171     b->atom1=d[0];
2172     b->atom2=d[1];
2173
2174     num_read++;
2175   }
2176
2177   return;
2178 }
2179 /*      END OF FUNCTION read_acceptors      */
2180
2181
2182 /************************************************************************/
2183 /*                  */
2184 /*      FUNCTION read_exclusions      */
2185 /*                  */
2186 /*   INPUTS:                */
2187 /*  fd - file descriptor for .psf file        */
2188 /*                  */
2189 /*  read_exclusions reads in the explicit non-bonded exclusions     */
2190 /*  from the .psf file.  This section is a little funky, so hang on.    */
2191 /*  Ok, first there is a list of atom indexes that is NumExclusions     */
2192 /*  long.  These are in some sense the atoms that will be exlcuded.     */
2193 /*  Following this list is a list of NumAtoms length that is a list     */
2194 /*  of indexes into the list of excluded atoms.  So an example.  Suppose*/
2195 /*  we have a 5 atom simulation with 3 explicit exclusions.  The .psf   */
2196 /*  file could look like:            */
2197 /*                  */
2198 /*  3!NNB                */
2199 /*  3 4 5                */
2200 /*  0 1 3 3 3              */
2201 /*                  */
2202 /*  This would mean that atom 1 has no explicit exclusions.  Atom 2     */
2203 /*  has an explicit exclusion with atom 3.  Atom 3 has an explicit      */
2204 /*  exclusion with atoms 4 AND 5.  And atoms 4 and 5 have no explicit   */
2205 /*  exclusions.  Got it!?!  I'm not sure who dreamed this up . . .      */
2206 /*                  */
2207 /************************************************************************/
2208
2209 void Molecule::read_exclusions(FILE *fd)
2210
2211 {
2212   int *exclusion_atoms;  //  Array of indexes of excluded atoms
2213   register int num_read=0;    //  Number fo exclusions read in
2214   int current_index;  //  Current index value
2215   int last_index;    //  the previous index value
2216   register int insert_index=0;  //  index of where we are in exlcusions array
2217
2218   /*  Allocate the array of exclusion structures and the array of */
2219   /*  exlcuded atom indexes          */
2220   exclusions      = new Exclusion[numExclusions];
2221   exclusion_atoms = new int[numExclusions];
2222
2223   if ( (exclusions == NULL) || (exclusion_atoms == NULL) )
2224   {
2225     NAMD_die("memory allocation failed in Molecule::read_exclusions");
2226   }
2227
2228   /*  First, read in the excluded atoms list      */
2229   for (num_read=0; num_read<numExclusions; num_read++)
2230   {
2231     /*  Read the atom number from the file. Subtract 1 to   */
2232     /*  convert the index from the 1 to NumAtoms used in the*/
2233     /*  file to the  0 to NumAtoms-1 that we need    */
2234     exclusion_atoms[num_read]=NAMD_read_int(fd, "IMPROPERS")-1;
2235
2236     /*  Check for an illegal index        */
2237     if (exclusion_atoms[num_read] >= numAtoms)
2238     {
2239       char err_msg[128];
2240
2241       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);
2242       NAMD_die(err_msg);
2243     }
2244   }
2245
2246   /*  Now, go through and read the list of NumAtoms pointers into */
2247   /*  the array that we just read in        */
2248   last_index=0;
2249
2250   for (num_read=0; num_read<numAtoms; num_read++)
2251   {
2252     /*  Read in the current index value      */
2253     current_index=NAMD_read_int(fd, "EXCLUSIONS");
2254
2255     /*  Check for an illegal pointer      */
2256     if (current_index>numExclusions)
2257     {
2258       char err_msg[128];
2259
2260       sprintf(err_msg, "EXCLUSION INDEX %d LARGER THAN NUMBER OF EXLCUSIONS %d IN PSF FILE, EXCLUSION #%d\n", 
2261          current_index+1, numExclusions, num_read);
2262       NAMD_die(err_msg);
2263     }
2264
2265     /*  Check to see if it matches the last index.  If so   */
2266     /*  than this atom has no exclusions.  If not, then     */
2267     /*  we have to build some exclusions      */
2268     if (current_index != last_index)
2269     {
2270       /*  This atom has some exclusions.  Loop from   */
2271       /*  the last_index to the current index.  This  */
2272       /*  will include how ever many exclusions this  */
2273       /*  atom has          */
2274       for (insert_index=last_index; 
2275            insert_index<current_index; insert_index++)
2276       {
2277         /*  Assign the two atoms involved.      */
2278         /*  The first one is our position in    */
2279         /*  the list, the second is based on    */
2280         /*  the pointer into the index list     */
2281         int a1 = num_read;
2282         int a2 = exclusion_atoms[insert_index];
2283         if ( a1 < a2 ) {
2284           exclusions[insert_index].atom1 = a1;
2285           exclusions[insert_index].atom2 = a2;
2286         } else if ( a2 < a1 ) {
2287           exclusions[insert_index].atom1 = a2;
2288           exclusions[insert_index].atom2 = a1;
2289         } else {
2290           char err_msg[128];
2291           sprintf(err_msg, "ATOM %d EXCLUDED FROM ITSELF IN PSF FILE\n", a1+1);
2292           NAMD_die(err_msg);
2293         }
2294       }
2295
2296       last_index=current_index;
2297     }
2298   }
2299
2300   /*  Free our temporary list of indexes        */
2301   delete [] exclusion_atoms;
2302
2303   return;
2304 }
2305 /*      END OF FUNCTION read_exclusions      */
2306
2307 /************************************************************************/
2308 /*      FUNCTION read_exclusions      */
2309 /*                  */
2310 /*   INPUTS:                */
2311 /*   int* atom_i - array of atom i indices    */
2312 /*   int* atom_j - array of atom j indices    */
2313 /*   int num_exclusion - length of array           */
2314 /*                  */
2315 /* JLai August 16th, 2012 */
2316 /************************************************************************/
2317 void Molecule::read_exclusions(int* atom_i, int* atom_j, int num_exclusion)
2318 {
2319     /*  Allocate the array of exclusion structures and the array of */
2320   /*  exlcuded atom indexes          */
2321   exclusions       = new Exclusion[num_exclusion];
2322   int loop_counter = 0;  
2323   int a=0;
2324   int b=0;
2325
2326   if ( (exclusions == NULL) )
2327   {
2328     NAMD_die("memory allocation failed in Molecule::read_exclusions");
2329   }
2330
2331   /* The following code only guarantees that exclusion.atom1 is < exclusion.atom2 */
2332   for (loop_counter = 0; loop_counter < num_exclusion; loop_counter++) {
2333         
2334         if ( (atom_i == NULL) || (atom_j == NULL) ) {
2335           NAMD_die("null pointer expection in Molecule::read_exclusions");
2336         }
2337
2338         a = atom_i[loop_counter];
2339         b = atom_j[loop_counter];
2340         if(a < b) {
2341                 exclusions[loop_counter].atom1 = a;
2342                 exclusions[loop_counter].atom2 = b;
2343         } else {
2344                 exclusions[loop_counter].atom1 = b;
2345                 exclusions[loop_counter].atom2 = a;
2346         }
2347         exclusionSet.add(Exclusion(exclusions[loop_counter].atom1,exclusions[loop_counter].atom2));
2348   }
2349
2350   if ( ! CkMyPe() ) {
2351     iout << iINFO << "ADDED " << num_exclusion << " EXPLICIT EXCLUSIONS: THIS VALUE WILL *NOT* BE ADDED TO THE STRUCTURE SUMMARY\n" << endi;
2352   }
2353
2354    return;
2355 }
2356 /*      END OF FUNCTION read_exclusions      */
2357
2358 /************************************************************************
2359  *
2360  *        FUNCTION read_lphosts
2361  *
2362  *   INPUTS:
2363  *  fd - file pointer to the .psf file
2364  *
2365  *  This function reads in the lone pair host section of the .psf file.
2366  *  Each lone pair is supported by 2 or 3 host atoms.
2367  *
2368  *  All lonepair specifications in a PSF are expected to have three
2369  *  associated values.  However, the meaning of these values depends
2370  *  on the lonepair type.  Nonetheless, for simplicity with the old
2371  *  code, the struct values are still called "distance", "angle",
2372  *  and "dihedral", even when this is not what the value signifies.
2373  *
2374  ************************************************************************/
2375 void Molecule::read_lphosts(FILE *fd)
2376 {
2377   char buffer[512];  // Buffer for reading from file
2378   char weight[8];    // Weighting type identified by string --
2379                      // unused/unsupported outside of CHARMM RTF
2380                      // so we read it, but ignore it.
2381   int numhosts;  // Refers to the number of atoms that support the
2382                  // given lone pair, either 2 or 3.
2383   int index;  // 1-based index into the stream of numbers (8 per line)
2384               // that indicates the start of each LP host record.
2385   int next_index = 1;  // Forecast next index value as an error check.
2386   int i, read_count;
2387   Real value1, value2, value3; 
2388
2389   // called only if numLphosts > 0
2390   lphosts = new Lphost[numLphosts];
2391   if (lphosts == NULL) {
2392     NAMD_die("memory allocation failed in Molecule::read_lphosts");
2393   }
2394   for (i = 0;  i < numLphosts;  i++) {
2395     NAMD_read_line(fd, buffer);
2396     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
2397     read_count=sscanf(buffer, "%d %d %6s %f %f %f",
2398         &numhosts, &index, weight, &value1, &value2, &value3);
2399     // The "weight" specification in PSF remains hardcoded as "F"
2400     // (for false) inside NAMD.  Again, no extant force field uses
2401     // lonepairs with weighted placement, so this can really only be
2402     // specified inside an RTF, which NAMD never reads anyway.
2403     if (read_count != 6 || index != next_index ||
2404         strcmp(weight,"F") != 0 || numhosts < 2 || 3 < numhosts) {
2405       char err_msg[128];
2406       sprintf(err_msg, "BAD FORMAT FOR LPHOST LINE %d IN PSF FILE LINE\n"
2407           "LINE=%s\n", i+1, buffer);
2408       NAMD_die(err_msg);
2409     }
2410     lphosts[i].numhosts = numhosts; // currently must be 2 or 3
2411     next_index += numhosts + 1;  // add 1 to account for LP index
2412     if (numhosts == 2) {
2413       lphosts[i].distance = value1;
2414       lphosts[i].angle = value2;
2415       lphosts[i].dihedral = 0.0;  // ignore value3
2416     }
2417     else {  // numhosts == 3
2418       lphosts[i].distance = value1;
2419       lphosts[i].angle = value2 * (M_PI/180);     // convert to radians
2420       lphosts[i].dihedral = value3 * (M_PI/180);  // convert to radians
2421     }
2422   }
2423   for (i = 0;  i < numLphosts;  i++) {
2424     // Subtract 1 to get 0-based atom index
2425     lphosts[i].atom1 = NAMD_read_int(fd, "LPHOSTS")-1;
2426     lphosts[i].atom2 = NAMD_read_int(fd, "LPHOSTS")-1;
2427     lphosts[i].atom3 = NAMD_read_int(fd, "LPHOSTS")-1;
2428     // For numhosts==2, set unused atom4 to atom1
2429     lphosts[i].atom4 = ( lphosts[i].numhosts == 3 ?
2430         NAMD_read_int(fd, "LPHOSTS")-1 : lphosts[i].atom1 );
2431   }
2432 }
2433 /*      END OF FUNCTION read_lphosts    */
2434
2435 /************************************************************************/
2436 /*                  */
2437 /*        FUNCTION read_anisos     */
2438 /*                  */
2439 /*   INPUTS:                */
2440 /*  fd - file pointer to the .psf file        */
2441 /*                  */
2442 /*  this function reads in the anisotropic terms section of .psf file. */
2443 /*                  */
2444 void Molecule::read_anisos(FILE *fd)
2445 {
2446   char buffer[512];  // Buffer for reading from file
2447   int numhosts, index, i, read_count;
2448   Real k11, k22, k33;
2449
2450   anisos = new Aniso[numAnisos];
2451   if (anisos == NULL)
2452   {
2453     NAMD_die("memory allocation failed in Molecule::read_anisos");
2454   }
2455   for (i = 0;  i < numAnisos;  i++)
2456   {
2457     NAMD_read_line(fd, buffer);
2458     if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') ) continue;
2459     read_count=sscanf(buffer, "%f %f %f", &k11, &k22, &k33);
2460     if (read_count != 3)
2461     {
2462       char err_msg[128];
2463       sprintf(err_msg, "BAD FORMAT FOR ANISO LINE %d IN PSF FILE LINE\n"
2464           "LINE=%s\n", i+1, buffer);
2465       NAMD_die(err_msg);
2466     }
2467     anisos[i].k11 = k11;
2468     anisos[i].k22 = k22;
2469     anisos[i].k33 = k33;
2470   }
2471   for (i = 0;  i < numAnisos;  i++) {
2472     anisos[i].atom1 = NAMD_read_int(fd, "ANISOS")-1;
2473     anisos[i].atom2 = NAMD_read_int(fd, "ANISOS")-1;
2474     anisos[i].atom3 = NAMD_read_int(fd, "ANISOS")-1;
2475     anisos[i].atom4 = NAMD_read_int(fd, "ANISOS")-1;
2476   }
2477 }
2478 /*      END OF FUNCTION read_anisos     */
2479
2480 //LCPO
2481 inline int getLCPOTypeAmber(char atomType[11], int numBonds) {
2482
2483   //Hydrogen
2484   if (atomType[0] == 'H' || atomType[0] == 'h') {
2485     return 0;
2486
2487   //Carbon
2488   } else if (atomType[0] == 'C' || atomType[0] == 'c') {
2489     if (//Sp3 Carbon
2490       //atomType[1] == 'T')// ||
2491       strcmp(atomType, "CT" )==0 )
2492       //strcmp(atomType, "CP1" )==0 ||
2493       //strcmp(atomType, "CP2" )==0 ||
2494       //strcmp(atomType, "CP3" )==0 ||
2495       //strcmp(atomType, "CS"  )==0 )
2496       {
2497       if (numBonds == 1)
2498         return 1;
2499       else if (numBonds == 2)
2500         return 2;
2501       else if (numBonds == 3)
2502         return 3;
2503       else if (numBonds == 4)
2504         return 4;
2505       else
2506         return 1;
2507
2508     } else {//Sp2 or other
2509       if (numBonds == 2)
2510         return 5;
2511       else if (numBonds == 3)
2512         return 6;
2513       else
2514         return 1;
2515     }
2516
2517   //Nitrogen
2518   } else if (atomType[0] == 'N' || atomType[0] == 'n') {
2519     if ( strcmp(atomType, "N3"  ) == 0 ) { //Sp3 Nitrogen
2520       if (numBonds == 1)
2521         return 11;
2522       else if (numBonds == 2)
2523         return 12;
2524       else if (numBonds == 3)
2525         return 13;
2526       else
2527         return 11;
2528
2529     } else {//SP2 Nitrogen
2530       if (numBonds == 1)
2531         return 14;
2532       else if (numBonds == 2)
2533         return 15;
2534       else if (numBonds == 3)
2535         return 16;
2536       else
2537         return 11;
2538     }
2539
2540   //Oxygen
2541   } else if (atomType[0] == 'O' || atomType[0] == 'o') {
2542
2543     if ( strcmp(atomType, "O" )==0) {//Sp2 Oxygen
2544       return 9;
2545     } else if (strcmp(atomType, "O2" )==0) {//Carboxylate Oxygen
2546       return 10;
2547     } else { // Sp3 Oxygen
2548       if (numBonds == 1)
2549         return 7;
2550       else if (numBonds == 2)
2551         return 8;
2552       else
2553         return 7;
2554     }
2555
2556   //Sulfur
2557   } else if (atomType[0] == 'S' || atomType[0] == 's') {
2558     if ( strcmp(atomType, "SH" )==0) { //Sulfur 1 neighbor
2559       return 17;
2560     } else {
2561       return 18;
2562     }
2563
2564   //Phosphorus
2565   } else if (atomType[0] == 'P' || atomType[0] == 'p') {
2566       if (numBonds == 3)
2567         return 19;
2568       else if (numBonds == 4)
2569         return 20;
2570       else
2571         return 19;
2572   } else if (atomType[0] == 'Z') { // ? just to agree with Amber mdread.f
2573     return 0;
2574   } else  if ( strcmp(atomType, "MG" )==0) { //Mg
2575     return 22;
2576   } else { // unknown atom type
2577     return 5;
2578   }
2579   return 5;
2580 } // getLCPOTypeAmber
2581
2582 inline int getLCPOTypeCharmm(char atomType[11], int numBonds) {
2583
2584   //Hydrogen
2585   if (atomType[0] == 'H') {
2586     return 0;
2587
2588   //Carbon
2589   } else if (atomType[0] == 'C') {
2590     if (//Sp3 Carbon
2591       atomType[1] == 'T' ||
2592       strcmp(atomType, "CP1" )==0 ||
2593       strcmp(atomType, "CP2" )==0 ||
2594       strcmp(atomType, "CP3" )==0 ||
2595       strcmp(atomType, "CS"  )==0 ) {
2596       if (numBonds == 1)
2597         return 1;
2598       else if (numBonds == 2)
2599         return 2;
2600       else if (numBonds == 3)
2601         return 3;
2602       else if (numBonds == 4)
2603         return 4;
2604       else
2605         return 1;
2606
2607     } else if (//Sp2
2608       strcmp(atomType, "C"   )==0 ||
2609       strcmp(atomType, "CA"  )==0 ||
2610       strcmp(atomType, "CC"  )==0 ||
2611       strcmp(atomType, "CD"  )==0 ||
2612       strcmp(atomType, "CN"  )==0 ||
2613       strcmp(atomType, "CY"  )==0 ||
2614       strcmp(atomType, "C3"  )==0 ||
2615       strcmp(atomType, "CE1" )==0 ||
2616       strcmp(atomType, "CE2" )==0 ||
2617       strcmp(atomType, "CST" )==0 ||
2618       strcmp(atomType, "CAP" )==0 ||
2619       strcmp(atomType, "COA" )==0 ||
2620       strcmp(atomType, "CPT" )==0 ||
2621       strcmp(atomType, "CPH1")==0 ||
2622       strcmp(atomType, "CPH2")==0
2623       ) {
2624       if (numBonds == 2)
2625         return 5;
2626       else if (numBonds == 3)
2627         return 6;
2628       else
2629         return 1;
2630     } else { // other Carbon
2631         return 1;
2632     }
2633
2634   //Nitrogen
2635   } else if (atomType[0] == 'N') {
2636     if (//Sp3 Nitrogen
2637       //strcmp(atomType, "N"   )==0 ||
2638       //strcmp(atomType, "NH1" )==0 ||
2639       //strcmp(atomType, "NH2" )==0 ||
2640       strcmp(atomType, "NH3" )==0 ||
2641       //strcmp(atomType, "NC2" )==0 ||
2642       //strcmp(atomType, "NY"  )==0 ||
2643       strcmp(atomType, "NP"  )==0
2644       ) {
2645       if (numBonds == 1)
2646         return 11;
2647       else if (numBonds == 2)
2648         return 12;
2649       else if (numBonds == 3)
2650         return 13;
2651       else
2652         return 11;
2653
2654     } else if (//SP2 Nitrogen
2655       strcmp(atomType, "NY"  )==0 || //
2656       strcmp(atomType, "NC2" )==0 || //
2657       strcmp(atomType, "N"   )==0 || //
2658       strcmp(atomType, "NH1" )==0 || //
2659       strcmp(atomType, "NH2" )==0 || //
2660       strcmp(atomType, "NR1" )==0 ||
2661       strcmp(atomType, "NR2" )==0 ||
2662       strcmp(atomType, "NR3" )==0 ||
2663       strcmp(atomType, "NPH" )==0 ||
2664       strcmp(atomType, "NC"  )==0
2665       ) {
2666       if (numBonds == 1)
2667         return 14;
2668       else if (numBonds == 2)
2669         return 15;
2670       else if (numBonds == 3)
2671         return 16;
2672       else
2673         return 11;
2674     } else { // other Nitrogen
2675       return 11;
2676     }
2677
2678   //Oxygen
2679   } else if (atomType[0] == 'O') {
2680     if (//Sp3 Oxygen
2681       strcmp(atomType, "OH1" )==0 ||
2682       strcmp(atomType, "OS"  )==0 ||
2683       strcmp(atomType, "OC"  )==0 || //
2684       strcmp(atomType, "OT"  )==0
2685       ) {
2686       if (numBonds == 1)
2687         return 7;
2688       else if (numBonds == 2)
2689         return 8;
2690       else
2691         return 7;
2692     } else if ( // Sp2 Oxygen
2693       strcmp(atomType, "O"   )==0 ||
2694       strcmp(atomType, "OB"  )==0 ||
2695       strcmp(atomType, "OST" )==0 ||
2696       strcmp(atomType, "OCA" )==0 ||
2697       strcmp(atomType, "OM"  )==0
2698       ) {
2699       return 9;
2700     } else if ( // SP1 Oxygen
2701       strcmp(atomType, "OC"  )==0
2702       ) {
2703       return 10;
2704     } else { // other Oxygen
2705       return 7;
2706     }
2707
2708   //Sulfur
2709   } else if (atomType[0] == 'S') {
2710       if (numBonds == 1)
2711         return 17;
2712       else
2713         return 18;
2714
2715   //Phosphorus
2716   } else if (atomType[0] == 'P') {
2717       if (numBonds == 3)
2718         return 19;
2719       else if (numBonds == 4)
2720         return 20;
2721       else
2722         return 19;
2723   } else { // unknown atom type
2724     return 5;
2725   }
2726   return 5;
2727 } // getLCPOTypeCharmm
2728
2729 //input type is Charmm/Amber/other
2730 //0 - Charmm/Xplor
2731 //1 - Amber
2732 //2 - Plugin
2733 //3 - Gromacs
2734 void Molecule::assignLCPOTypes(int inputType) {
2735   int *heavyBonds = new int[numAtoms];
2736   for (int i = 0; i < numAtoms; i++)
2737     heavyBonds[i] = 0;
2738   for (int i = 0; i < numBonds; i++ ) {
2739     Bond curBond = bonds[i];
2740     int a1 = bonds[i].atom1;
2741     int a2 = bonds[i].atom2;
2742     if (atoms[a1].mass > 2.f && atoms[a2].mass > 2.f) {
2743       heavyBonds[a1]++;
2744       heavyBonds[a2]++;
2745     }
2746   }
2747
2748   lcpoParamType = new int[numAtoms];
2749
2750   int warning = 0;
2751   for (int i = 0; i < numAtoms; i++) {
2752     //print vdw_type and numbonds
2753
2754     if (inputType == 1) { // Amber
2755       lcpoParamType[i] = getLCPOTypeAmber(atomNames[i].atomtype, heavyBonds[i]);
2756     } else { // Charmm
2757       lcpoParamType[i] = getLCPOTypeCharmm(atomNames[i].atomtype, heavyBonds[i]);
2758     }
2759 /*
2760     CkPrintf("%d MOL: ATOM[%05d] = { %4s %d } : %d\n",
2761       inputType,
2762       i+1,
2763       atomNames[i].atomtype,
2764       heavyBonds[i],
2765       lcpoParamType[i]
2766       );
2767 */
2768     if ( atoms[i].mass < 1.5 && lcpoParamType[i] != 0 ) {
2769       if (atoms[i].status & LonepairAtom) {
2770         warning |= LonepairAtom;
2771         lcpoParamType[i] = 0;  // reset LCPO type for LP
2772       }
2773       else if (atoms[i].status & DrudeAtom) {
2774         warning |= DrudeAtom;
2775         lcpoParamType[i] = 0;  // reset LCPO type for Drude
2776       }
2777       else {
2778         CkPrintf("ERROR in Molecule::assignLCPOTypes(): "
2779             "Light atom given heavy atom LCPO type.\n");
2780       }
2781     }
2782
2783     //CkPrintf("VDW_TYPE %02d %4s\n", atoms[i].vdw_type, atomNames[i].atomtype);
2784   } // for atoms
2785
2786   if (warning & LonepairAtom) {
2787     iout << iWARN << "LONE PAIRS TO BE IGNORED BY SASA\n" << endi;
2788   }
2789   if (warning & DrudeAtom) {
2790     iout << iWARN << "DRUDE PARTICLES TO BE IGNORED BY SASA\n" << endi;
2791   }
2792
2793   delete [] heavyBonds;
2794
2795 } // buildLCPOTable
2796
2797 void Molecule::plgLoadAtomBasics(molfile_atom_t *atomarray){
2798     atoms = new Atom[numAtoms];
2799     atomNames = new AtomNameInfo[numAtoms];
2800     if(simParams->genCompressedPsf) {
2801         atomSegResids = new AtomSegResInfo[numAtoms];
2802     }    
2803     hydrogenGroup.resize(0);
2804
2805     ResidueLookupElem *tmpResLookup = resLookup;
2806
2807     for(int i=0; i<numAtoms; i++) {
2808         int reslength = strlen(atomarray[i].resname)+1;
2809         int namelength = strlen(atomarray[i].name)+1;
2810         int typelength = strlen(atomarray[i].type)+1;
2811         atomNames[i].resname = nameArena->getNewArray(reslength);
2812         atomNames[i].atomname = nameArena->getNewArray(namelength);
2813         atomNames[i].atomtype = nameArena->getNewArray(typelength);
2814         strcpy(atomNames[i].resname, atomarray[i].resname);
2815         strcpy(atomNames[i].atomname, atomarray[i].name);
2816         strcpy(atomNames[i].atomtype, atomarray[i].type);
2817
2818         atoms[i].mass = atomarray[i].mass;
2819         atoms[i].charge = atomarray[i].charge;
2820         atoms[i].status = UnknownAtom;
2821
2822         //add this atom to residue lookup table
2823         if(tmpResLookup) {
2824             tmpResLookup = tmpResLookup->append(atomarray[i].segid, atomarray[i].resid, i);
2825         }
2826
2827         if(atomSegResids) { //for compressing molecule information
2828             AtomSegResInfo *one = atomSegResids + i;
2829             memcpy(one->segname, atomarray[i].segid, strlen(atomarray[i].segid)+1);
2830             one->resid = atomarray[i].resid;
2831         }
2832         //Determine the type of the atom
2833         if ( simParams->ignoreMass ) {
2834         }else if(atoms[i].mass <= 0.05) {
2835             atoms[i].status |= LonepairAtom;
2836         }else if(atoms[i].mass < 1.0) {
2837             atoms[i].status |= DrudeAtom;
2838         }else if(atoms[i].mass <= 3.5) {
2839             atoms[i].status |= HydrogenAtom;
2840         }else if((atomNames[i].atomname[0] == 'O') &&
2841                  (atoms[i].mass>=14.0) && (atoms[i].mass<=18.0)){
2842             atoms[i].status |= OxygenAtom;
2843         }
2844         //Look up the vdw constants for this atom
2845         params->assign_vdw_index(atomNames[i].atomtype, &atoms[i]);
2846     }
2847 }
2848
2849 void Molecule::plgLoadBonds(int *from, int *to){
2850     bonds = new Bond[numBonds];
2851     char atom1name[11];
2852     char atom2name[11];
2853     int realNumBonds = 0;
2854     for(int i=0; i<numBonds; i++) {
2855         Bond *thisBond = bonds+realNumBonds;
2856         thisBond->atom1 = from[i]-1;
2857         thisBond->atom2 = to[i]-1;
2858         /* Get the atom type for the two atoms.
2859          * When we query the parameter object, we
2860          * need to send the atom type that is alphabetically
2861          * first as atom 1.
2862          */
2863         if(strcasecmp(atomNames[thisBond->atom1].atomtype,
2864                       atomNames[thisBond->atom2].atomtype)<0) {
2865             strcpy(atom1name, atomNames[thisBond->atom1].atomtype);
2866             strcpy(atom2name, atomNames[thisBond->atom2].atomtype);
2867         }else{
2868             strcpy(atom2name, atomNames[thisBond->atom1].atomtype);
2869             strcpy(atom1name, atomNames[thisBond->atom2].atomtype);
2870         }
2871         params->assign_bond_index(atom1name, atom2name, thisBond);
2872
2873         //Make sure this isn't a fake bond meant for shake in x-plor
2874         Real k, x0;
2875         params->get_bond_params(&k, &x0, thisBond->bond_type);
2876         if (is_lonepairs_psf) {
2877             //need to retain Lonepair bonds for Drude
2878             if(k!=0. || is_lp(thisBond->atom1) || 
2879                is_lp(thisBond->atom2)) {               
2880                 realNumBonds++;
2881             }
2882         }else{
2883             if(k != 0.) realNumBonds++;
2884         }
2885     }
2886
2887     if(numBonds != realNumBonds) {
2888         iout << iWARN << "Ignored" << numBonds-realNumBonds <<
2889             "bonds with zero force constants.\n" <<endi;
2890         iout << iWARN << "Will get H-H distance in rigid H20 from H-O-H angle.\n" <<endi;
2891     }
2892     numBonds = realNumBonds;
2893 }
2894
2895 void Molecule::plgLoadAngles(int *plgAngles)
2896 {    
2897     char atom1name[11];
2898     char atom2name[11];
2899     char atom3name[11];
2900
2901     angles=new Angle[numAngles];
2902     int *atomid = plgAngles;
2903     int numRealAngles = 0;
2904     for(int i=0; i<numAngles; i++) {
2905         Angle *thisAngle = angles+numRealAngles;
2906         thisAngle->atom1 = atomid[0]-1;
2907         thisAngle->atom2 = atomid[1]-1;
2908         thisAngle->atom3 = atomid[2]-1;
2909         atomid += 3;
2910
2911         if(strcasecmp(atomNames[thisAngle->atom1].atomtype,
2912                       atomNames[thisAngle->atom2].atomtype)<0) {
2913             strcpy(atom1name, atomNames[thisAngle->atom1].atomtype);
2914             strcpy(atom2name, atomNames[thisAngle->atom2].atomtype);
2915             strcpy(atom3name, atomNames[thisAngle->atom3].atomtype);
2916         }else{
2917             strcpy(atom1name, atomNames[thisAngle->atom3].atomtype);
2918             strcpy(atom2name, atomNames[thisAngle->atom2].atomtype);
2919             strcpy(atom3name, atomNames[thisAngle->atom1].atomtype);
2920         }
2921
2922         params->assign_angle_index(atom1name, atom2name, atom3name,
2923                                 thisAngle, simParams->alchOn ? -1 : 0);
2924         if ( thisAngle->angle_type == -1 ) {
2925           iout << iWARN << "ALCHEMY MODULE WILL REMOVE ANGLE OR RAISE ERROR\n"
2926                << endi;
2927         }
2928
2929         Real k, t0, k_ub, r_ub;
2930         if ( thisAngle->angle_type == -1 ) { k = -1.;  k_ub = -1.; } else
2931         params->get_angle_params(&k, &t0, &k_ub, &r_ub, thisAngle->angle_type);
2932         if(k!=0. || k_ub!=0.) numRealAngles++;
2933     }
2934
2935     if(numAngles != numRealAngles) {
2936         iout << iWARN << "Ignored" << numAngles-numRealAngles << 
2937             " angles with zero force constants.\n" << endi; 
2938     }
2939     numAngles = numRealAngles;
2940 }
2941
2942 void Molecule::plgLoadDihedrals(int *plgDihedrals)
2943 {
2944     std::map< std::string, int > cache;
2945
2946     int lastAtomIds[4];
2947     int multiplicity = 1; //multiplicity of the current bond
2948
2949     lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
2950     dihedrals = new Dihedral[numDihedrals];
2951     int numRealDihedrals = 0;
2952     int *atomid = plgDihedrals;
2953     for(int i=0; i<numDihedrals; i++, atomid+=4) {
2954         Dihedral *thisDihedral = dihedrals + numRealDihedrals;
2955         Bool duplicate_bond = TRUE;
2956         for(int j=0; j<4; j++) {
2957             if(atomid[j] != lastAtomIds[j]) {
2958                 duplicate_bond = FALSE;
2959             }
2960             lastAtomIds[j] = atomid[j];
2961         }
2962
2963         if(duplicate_bond) {
2964             multiplicity++;
2965             if(multiplicity==2) {
2966                 numMultipleDihedrals++;
2967             }
2968         }else{
2969             multiplicity=1;
2970             numRealDihedrals++;
2971         }
2972
2973         thisDihedral->atom1 = atomid[0]-1;
2974         thisDihedral->atom2 = atomid[1]-1;
2975         thisDihedral->atom3 = atomid[2]-1;
2976         thisDihedral->atom4 = atomid[3]-1;
2977
2978       char query[128];
2979       sprintf(query,"%10s :: %10s :: %10s :: %10s :: %d",
2980         atomNames[atomid[0]-1].atomtype,
2981         atomNames[atomid[1]-1].atomtype,
2982         atomNames[atomid[2]-1].atomtype,
2983         atomNames[atomid[3]-1].atomtype,
2984         multiplicity);
2985       auto search = cache.find(query);
2986       if ( search != cache.end() ) { 
2987         thisDihedral->dihedral_type = search->second;
2988       } else {
2989         char atom1name[11];
2990         char atom2name[11];
2991         char atom3name[11];
2992         char atom4name[11];
2993         strcpy(atom1name, atomNames[atomid[0]-1].atomtype);
2994         strcpy(atom2name, atomNames[atomid[1]-1].atomtype);
2995         strcpy(atom3name, atomNames[atomid[2]-1].atomtype);
2996         strcpy(atom4name, atomNames[atomid[3]-1].atomtype);
2997
2998         params->assign_dihedral_index(atom1name, atom2name,
2999                                       atom3name, atom4name, thisDihedral,
3000                                       multiplicity, simParams->alchOn ? -1 : 0);
3001         if ( thisDihedral->dihedral_type == -1 ) {
3002           iout << iWARN << "ALCHEMY MODULE WILL REMOVE DIHEDRAL OR RAISE ERROR\n"
3003                << endi;
3004         }
3005         cache[query] = thisDihedral->dihedral_type;
3006       }
3007     }
3008
3009     numDihedrals = numRealDihedrals;
3010 }
3011
3012 void Molecule::plgLoadImpropers(int *plgImpropers)
3013 {
3014     char atom1name[11];
3015     char atom2name[11];
3016     char atom3name[11];
3017     char atom4name[11];
3018     int lastAtomIds[4];
3019     int multiplicity = 1; //multiplicity of the current bond
3020
3021     lastAtomIds[0]=lastAtomIds[1]=lastAtomIds[2]=lastAtomIds[3]=-1;
3022     impropers = new Improper[numImpropers];
3023     int numRealImpropers = 0;
3024     int *atomid = plgImpropers;
3025     for(int i=0; i<numImpropers; i++, atomid+=4) {
3026         Improper *thisImproper = impropers + numRealImpropers;
3027         Bool duplicate_bond = TRUE;
3028         for(int j=0; j<4; j++) {
3029             if(atomid[j] != lastAtomIds[j]) {
3030                 duplicate_bond = FALSE;
3031             }
3032             lastAtomIds[j] = atomid[j];
3033         }
3034
3035         strcpy(atom1name, atomNames[atomid[0]-1].atomtype);
3036         strcpy(atom2name, atomNames[atomid[1]-1].atomtype);
3037         strcpy(atom3name, atomNames[atomid[2]-1].atomtype);
3038         strcpy(atom4name, atomNames[atomid[3]-1].atomtype);
3039
3040         if(duplicate_bond) {
3041             multiplicity++;
3042             if(multiplicity==2) {
3043                 numMultipleImpropers++;
3044             }
3045         }else{
3046             multiplicity=1;
3047             numRealImpropers++;
3048         }
3049
3050         thisImproper->atom1 = atomid[0]-1;
3051         thisImproper->atom2 = atomid[1]-1;
3052         thisImproper->atom3 = atomid[2]-1;
3053         thisImproper->atom4 = atomid[3]-1;
3054
3055         params->assign_improper_index(atom1name, atom2name,
3056                                       atom3name, atom4name, thisImproper,
3057                                       multiplicity);
3058     }
3059
3060     numImpropers = numRealImpropers;
3061 }
3062
3063 void Molecule::plgLoadCrossterms(int *plgCterms)
3064 {
3065     char atom1name[11];
3066     char atom2name[11];
3067     char atom3name[11];
3068     char atom4name[11];
3069     char atom5name[11];
3070     char atom6name[11];
3071     char atom7name[11];
3072     char atom8name[11];
3073     int lastAtomIds[8];    
3074
3075     for(int i=0; i<8; i++)
3076         lastAtomIds[i]=-1;
3077     
3078     crossterms = new Crossterm[numCrossterms];
3079     int numRealCrossterms = 0;
3080     int *atomid = plgCterms;
3081     for(int i=0; i<numCrossterms; i++, atomid+=8) {
3082         Crossterm *thisCrossterm = crossterms + numRealCrossterms;
3083         Bool duplicate_bond = TRUE;
3084         for(int j=0; j<8; j++) {
3085             if(atomid[j] != lastAtomIds[j]) {
3086                 duplicate_bond = FALSE;
3087             }
3088             lastAtomIds[j] = atomid[j];
3089         }
3090
3091         strcpy(atom1name, atomNames[atomid[0]-1].atomtype);
3092         strcpy(atom2name, atomNames[atomid[1]-1].atomtype);
3093         strcpy(atom3name, atomNames[atomid[2]-1].atomtype);
3094         strcpy(atom4name, atomNames[atomid[3]-1].atomtype);
3095         strcpy(atom5name, atomNames[atomid[4]-1].atomtype);
3096         strcpy(atom6name, atomNames[atomid[5]-1].atomtype);
3097         strcpy(atom7name, atomNames[atomid[6]-1].atomtype);
3098         strcpy(atom8name, atomNames[atomid[7]-1].atomtype);
3099
3100         if(duplicate_bond) {
3101             iout << iWARN <<"Duplicate cross-term detected.\n" << endi;
3102         } else
3103             numRealCrossterms++;
3104
3105         thisCrossterm->atom1 = atomid[0]-1;
3106         thisCrossterm->atom2 = atomid[1]-1;
3107         thisCrossterm->atom3 = atomid[2]-1;
3108         thisCrossterm->atom4 = atomid[3]-1;
3109         thisCrossterm->atom5 = atomid[4]-1;
3110         thisCrossterm->atom6 = atomid[5]-1;
3111         thisCrossterm->atom7 = atomid[6]-1;
3112         thisCrossterm->atom8 = atomid[7]-1;
3113
3114         params->assign_crossterm_index(atom1name, atom2name,
3115                                        atom3name, atom4name, atom5name,
3116                                        atom6name, atom7name, atom8name,
3117                                        thisCrossterm);
3118     }
3119
3120     numCrossterms = numRealCrossterms;
3121 }
3122
3123 void Molecule::setOccupancyData(molfile_atom_t *atomarray){
3124     occupancy = new float[numAtoms];
3125     for(int i=0; i<numAtoms; i++) {
3126         occupancy[i] = atomarray[i].occupancy;
3127     }
3128 }
3129
3130 void Molecule::setBFactorData(molfile_atom_t *atomarray){
3131     bfactor = new float[numAtoms];
3132     for(int i=0; i<numAtoms; i++) {
3133         bfactor[i] = atomarray[i].bfactor;
3134     }
3135 }
3136
3137     /************************************************************************/
3138     /*                  */
3139     /*      FUNCTION build_lists_by_atom      */
3140     /*                  */
3141     /*  This function builds O(NumAtoms) arrays that store the bonds,   */
3142     /*  angles, dihedrals, and impropers, that each atom is involved in.    */
3143     /*  This is a space hog, but VERY fast.  This will certainly have to    */
3144     /*  change to make things scalable in memory, but for now, speed is the */
3145     /*  thing!                */
3146     /*                  */
3147     /************************************************************************/
3148     void Molecule::build_lists_by_atom()       
3149     {
3150        register int i;      //  Loop counter
3151
3152        register int numFixedAtoms = this->numFixedAtoms;
3153        // if we want forces on fixed atoms then just pretend
3154        // there are none for the purposes of this routine;
3155        if ( simParams->fixedAtomsForces ) numFixedAtoms = 0;
3156
3157 //fepb
3158 //     int numFepInitial = this->numFepInitial;
3159 //     int numFepFinal = this->numFepFinal;
3160 //fepe
3161        tmpArena = new ObjectArena<int32>;
3162        bondsWithAtom = new int32 *[numAtoms];
3163        cluster = new int32 [numAtoms];
3164        clusterSize = new int32 [numAtoms];
3165
3166        bondsByAtom = new int32 *[numAtoms];
3167        anglesByAtom = new int32 *[numAtoms];
3168        dihedralsByAtom = new int32 *[numAtoms];
3169        impropersByAtom = new int32 *[numAtoms];
3170        crosstermsByAtom = new int32 *[numAtoms];
3171
3172        exclusionsByAtom = new int32 *[numAtoms];
3173        fullExclusionsByAtom = new int32 *[numAtoms];
3174        modExclusionsByAtom = new int32 *[numAtoms];
3175
3176        // JLai
3177        gromacsPairByAtom = new int32 *[numAtoms];
3178        // End of JLai
3179
3180        int32 *byAtomSize = new int32[numAtoms];
3181
3182        const int pair_self = 
3183          simParams->pairInteractionOn ? simParams->pairInteractionSelf : 0;
3184
3185        DebugM(3,"Building bond lists.\n");
3186     
3187        //  Build the bond lists
3188        for (i=0; i<numAtoms; i++)
3189        {
3190          byAtomSize[i] = 0;
3191        }
3192        for (i=0; i<numRealBonds; i++)
3193        {
3194          byAtomSize[bonds[i].atom1]++;
3195          byAtomSize[bonds[i].atom2]++;
3196        }
3197        for (i=0; i<numAtoms; i++)
3198        {
3199          bondsWithAtom[i] = tmpArena->getNewArray(byAtomSize[i]+1);
3200          bondsWithAtom[i][byAtomSize[i]] = -1;
3201          byAtomSize[i] = 0;
3202        }
3203        for (i=0; i<numRealBonds; i++)
3204        {
3205          int a1 = bonds[i].atom1;
3206          int a2 = bonds[i].atom2;
3207          bondsWithAtom[a1][byAtomSize[a1]++] = i;
3208          bondsWithAtom[a2][byAtomSize[a2]++] = i;
3209        }
3210
3211         
3212        // Updates all bond, angle, dihedral, improper and crossterm
3213        // to reflect the QM region (which can't have any of there terms)
3214        if (simParams->qmForcesOn) {
3215            
3216            DebugM(3,"Calculating exclusions for QM simulation.\n");
3217            build_exclusions();
3218            
3219            delete_qm_bonded() ;
3220            
3221            DebugM(3,"Re-Building bond lists.\n");
3222            
3223            // We re-calculate the bondsWithAtom list for cluster 
3224            // info calculation below.
3225            for (i=0; i<numAtoms; i++)
3226            {
3227              byAtomSize[i] = 0;
3228            }
3229            for (i=0; i<numRealBonds; i++)
3230            {
3231              byAtomSize[bonds[i].atom1]++;
3232              byAtomSize[bonds[i].atom2]++;
3233            }
3234            for (i=0; i<numAtoms; i++)
3235            {
3236              bondsWithAtom[i][byAtomSize[i]] = -1;
3237              byAtomSize[i] = 0;
3238            }
3239            for (i=0; i<numRealBonds; i++)
3240            {
3241              int a1 = bonds[i].atom1;
3242              int a2 = bonds[i].atom2;
3243              bondsWithAtom[a1][byAtomSize[a1]++] = i;
3244              bondsWithAtom[a2][byAtomSize[a2]++] = i;
3245            }
3246        }
3247         
3248        //  Build cluster information (contiguous clusters)
3249        for (i=0; i<numAtoms; i++) {
3250          cluster[i] = i;
3251        }
3252        for (i=0; i<numAtoms; i++) {
3253          int ci = i;
3254          while ( cluster[ci] != ci ) ci = cluster[ci];
3255          for ( int32 *b = bondsWithAtom[i]; *b != -1; ++b ) {
3256            int a = bonds[*b].atom1;
3257            if ( a == i ) a = bonds[*b].atom2;
3258            if ( a > i ) {
3259              int ca = a;
3260              while ( cluster[ca] != ca ) ca = cluster[ca];
3261              if ( ca > ci ) cluster[ca] = cluster[ci];
3262              else cluster[ci] = cluster[ca];
3263            }
3264          }
3265        }
3266        while ( 1 ) {
3267          int allok = 1;
3268          for (i=0; i<numAtoms; i++) {
3269            int ci = cluster[i];
3270            if ( cluster[ci] != ci ) {
3271              allok = 0;
3272              cluster[i] = cluster[ci];
3273            }
3274          }
3275          if ( allok ) break;
3276        }
3277        
3278        for (i=0; i<numAtoms; i++) {
3279          clusterSize[i] = 0;
3280        }       
3281        for (i=0; i<numAtoms; i++) {           
3282          clusterSize[cluster[i]] += 1;
3283        }
3284
3285 /*
3286        //Getting number of clusters for debugging
3287        int numClusters=0;
3288        for(int i=0; i<numAtoms; i++){
3289            if(clusterSize[i]!=0) numClusters++;
3290        }
3291        printf("Num of clusters: %d\n", numClusters);
3292 */
3293
3294        //  Build the bond lists
3295        for (i=0; i<numAtoms; i++)
3296        {
3297          byAtomSize[i] = 0;
3298        }
3299        numCalcBonds = 0;
3300        for (i=0; i<numBonds; i++)
3301        {
3302          if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3303                             && fixedAtomFlags[bonds[i].atom2] ) continue;
3304    
3305          if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
3306          byAtomSize[bonds[i].atom1]++;
3307          numCalcBonds++;
3308        }
3309        for (i=0; i<numAtoms; i++)
3310        {
3311          bondsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3312          bondsByAtom[i][byAtomSize[i]] = -1;
3313          byAtomSize[i] = 0;
3314        }
3315        for (i=0; i<numBonds; i++)
3316        {
3317          if ( numFixedAtoms && fixedAtomFlags[bonds[i].atom1]
3318                             && fixedAtomFlags[bonds[i].atom2] ) continue;
3319          if ( pair_self && fepAtomFlags[bonds[i].atom1] != 1) continue;
3320          int a1 = bonds[i].atom1;
3321          bondsByAtom[a1][byAtomSize[a1]++] = i;
3322        }
3323        for (i=0; i<numBonds; ++i) {
3324          int a1 = bonds[i].atom1;
3325          int a2 = bonds[i].atom2;
3326          int j;
3327          if ( a1 == a2 ) {
3328            char buff[512];
3329            sprintf(buff,"Atom %d is bonded to itself", a1+1);
3330            NAMD_die(buff);
3331          }
3332          for ( j = 0; j < byAtomSize[a1]; ++j ) {
3333            int b = bondsByAtom[a1][j];
3334            int ba1 = bonds[b].atom1;
3335            int ba2 = bonds[b].atom2;
3336            if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3337              char buff[512];
3338              sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3339              NAMD_die(buff);
3340            }
3341          }
3342          for ( j = 0; j < byAtomSize[a2]; ++j ) {
3343            int b = bondsByAtom[a2][j];
3344            int ba1 = bonds[b].atom1;
3345            int ba2 = bonds[b].atom2;
3346            if ( b != i && ( (ba1==a1 && ba2==a2) || (ba1==a2 && ba2==a1) ) ) {
3347              char buff[512];
3348              sprintf(buff,"Duplicate bond from atom %d to atom %d", a1+1, a2+1);
3349              NAMD_die(buff);
3350            }
3351          }
3352        }
3353        
3354        DebugM(3,"Building angle lists.\n");
3355     
3356        //  Build the angle lists
3357        for (i=0; i<numAtoms; i++)
3358        {
3359          byAtomSize[i] = 0;
3360        }
3361        numCalcAngles = 0;
3362        for (i=0; i<numAngles; i++)
3363        {
3364          if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3365                             && fixedAtomFlags[angles[i].atom2]
3366                             && fixedAtomFlags[angles[i].atom3] ) continue;
3367          if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
3368          byAtomSize[angles[i].atom1]++;
3369          numCalcAngles++;
3370        }
3371        for (i=0; i<numAtoms; i++)
3372        {
3373          anglesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3374          anglesByAtom[i][byAtomSize[i]] = -1;
3375          byAtomSize[i] = 0;
3376        }
3377        for (i=0; i<numAngles; i++)
3378        {
3379          if ( pair_self && fepAtomFlags[angles[i].atom1] != 1) continue;
3380          if ( numFixedAtoms && fixedAtomFlags[angles[i].atom1]
3381                             && fixedAtomFlags[angles[i].atom2]
3382                             && fixedAtomFlags[angles[i].atom3] ) continue;
3383          int a1 = angles[i].atom1;
3384          anglesByAtom[a1][byAtomSize[a1]++] = i;
3385        }
3386        
3387        DebugM(3,"Building improper lists.\n");
3388     
3389        //  Build the improper lists
3390        for (i=0; i<numAtoms; i++)
3391        {
3392          byAtomSize[i] = 0;
3393        }
3394        numCalcImpropers = 0;
3395        for (i=0; i<numImpropers; i++)
3396        {
3397          if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3398                             && fixedAtomFlags[impropers[i].atom2]
3399                             && fixedAtomFlags[impropers[i].atom3]
3400                             && fixedAtomFlags[impropers[i].atom4] ) continue;
3401          if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
3402          byAtomSize[impropers[i].atom1]++;
3403          numCalcImpropers++;
3404        }
3405        for (i=0; i<numAtoms; i++)
3406        {
3407          impropersByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3408          impropersByAtom[i][byAtomSize[i]] = -1;
3409          byAtomSize[i] = 0;
3410        }
3411        for (i=0; i<numImpropers; i++)
3412        {
3413          if ( numFixedAtoms && fixedAtomFlags[impropers[i].atom1]
3414                             && fixedAtomFlags[impropers[i].atom2]
3415                             && fixedAtomFlags[impropers[i].atom3]
3416                             && fixedAtomFlags[impropers[i].atom4] ) continue;
3417          if ( pair_self && fepAtomFlags[impropers[i].atom1] != 1) continue;
3418          int a1 = impropers[i].atom1;
3419          impropersByAtom[a1][byAtomSize[a1]++] = i;
3420        }
3421        
3422        DebugM(3,"Building dihedral lists.\n");
3423     
3424        //  Build the dihedral lists
3425        for (i=0; i<numAtoms; i++)
3426        {
3427          byAtomSize[i] = 0;
3428        }
3429        numCalcDihedrals = 0;
3430        for (i=0; i<numDihedrals; i++)
3431        {
3432          if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3433                             && fixedAtomFlags[dihedrals[i].atom2]
3434                             && fixedAtomFlags[dihedrals[i].atom3]
3435                             && fixedAtomFlags[dihedrals[i].atom4] ) continue;
3436          if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
3437          byAtomSize[dihedrals[i].atom1]++;
3438          numCalcDihedrals++;
3439        }
3440        for (i=0; i<numAtoms; i++)
3441        {
3442          dihedralsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3443          dihedralsByAtom[i][byAtomSize[i]] = -1;
3444          byAtomSize[i] = 0;
3445        }
3446        for (i=0; i<numDihedrals; i++)
3447        {
3448          if ( numFixedAtoms && fixedAtomFlags[dihedrals[i].atom1]
3449                             && fixedAtomFlags[dihedrals[i].atom2]
3450                             && fixedAtomFlags[dihedrals[i].atom3]
3451                             && fixedAtomFlags[dihedrals[i].atom4] ) continue;
3452          if ( pair_self && fepAtomFlags[dihedrals[i].atom1] != 1) continue;
3453          int a1 = dihedrals[i].atom1;
3454          dihedralsByAtom[a1][byAtomSize[a1]++] = i;
3455        }
3456     
3457        DebugM(3,"Building crossterm lists.\n");
3458     
3459        //  Build the crossterm lists
3460        for (i=0; i<numAtoms; i++)
3461        {
3462          byAtomSize[i] = 0;
3463        }
3464        numCalcCrossterms = 0;
3465        for (i=0; i<numCrossterms; i++)
3466        {
3467          if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3468                             && fixedAtomFlags[crossterms[i].atom2]
3469                             && fixedAtomFlags[crossterms[i].atom3]
3470                             && fixedAtomFlags[crossterms[i].atom4]
3471                             && fixedAtomFlags[crossterms[i].atom5]
3472                             && fixedAtomFlags[crossterms[i].atom6]
3473                             && fixedAtomFlags[crossterms[i].atom7]
3474                             && fixedAtomFlags[crossterms[i].atom8] ) continue;
3475          if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
3476          byAtomSize[crossterms[i].atom1]++;
3477          numCalcCrossterms++;
3478        }
3479        for (i=0; i<numAtoms; i++)
3480        {
3481          crosstermsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3482          crosstermsByAtom[i][byAtomSize[i]] = -1;
3483          byAtomSize[i] = 0;
3484        }
3485        for (i=0; i<numCrossterms; i++)
3486        {
3487          if ( numFixedAtoms && fixedAtomFlags[crossterms[i].atom1]
3488                             && fixedAtomFlags[crossterms[i].atom2]
3489                             && fixedAtomFlags[crossterms[i].atom3]
3490                             && fixedAtomFlags[crossterms[i].atom4]
3491                             && fixedAtomFlags[crossterms[i].atom5]
3492                             && fixedAtomFlags[crossterms[i].atom6]
3493                             && fixedAtomFlags[crossterms[i].atom7]
3494                             && fixedAtomFlags[crossterms[i].atom8] ) continue;
3495          if ( pair_self && fepAtomFlags[crossterms[i].atom1] != 1) continue;
3496          int a1 = crossterms[i].atom1;
3497          crosstermsByAtom[a1][byAtomSize[a1]++] = i;
3498        }
3499
3500        // DRUDE: init lphostIndexes array
3501        if (is_lonepairs_psf) {
3502          // allocate lone pair host index array only if we need it!
3503          DebugM(3,"Initializing lone pair host index array.\n");
3504          lphostIndexes = new int32[numAtoms];
3505          for (i = 0;  i < numAtoms;  i++) {
3506            lphostIndexes[i] = -1;
3507          }
3508          for (i = 0;  i < numLphosts;  i++) {
3509            int32 index = lphosts[i].atom1;
3510            lphostIndexes[index] = i;
3511          }
3512        }
3513        // DRUDE
3514     
3515        // JLai
3516        DebugM(3,"Building gromacsPair lists.\n");
3517     
3518        //  Build the gromacsPair lists
3519        for (i=0; i<numAtoms; i++)
3520        {
3521          byAtomSize[i] = 0;
3522        }
3523        numCalcLJPair = 0;
3524        for (i=0; i<numLJPair; i++)
3525        {
3526          if ( numFixedAtoms && fixedAtomFlags[gromacsPair[i].atom1]
3527                             && fixedAtomFlags[gromacsPair[i].atom2] ) continue;
3528          if ( pair_self && fepAtomFlags[gromacsPair[i].atom1] != 1) continue;
3529          byAtomSize[gromacsPair[i].atom1]++;
3530          numCalcLJPair++;
3531        }
3532        for (i=0; i<numAtoms; i++)
3533        {
3534          gromacsPairByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3535          gromacsPairByAtom[i][byAtomSize[i]] = -1;
3536          byAtomSize[i] = 0;
3537        }
3538        for (i=0; i<numLJPair; i++)
3539        {
3540            if ( numFixedAtoms && fixedAtomFlags[gromacsPair[i].atom1]
3541                 && fixedAtomFlags[gromacsPair[i].atom2] ) continue;
3542            if ( pair_self && fepAtomFlags[gromacsPair[i].atom1] != 1) continue;
3543            int a1 = gromacsPair[i].atom1;
3544            gromacsPairByAtom[a1][byAtomSize[a1]++] = i;
3545            }
3546
3547        // End of JLai
3548
3549        DebugM(3,"Building exclusion data.\n");
3550     
3551        //  Build the arrays of exclusions for each atom
3552        if (! simParams->qmForcesOn)
3553        build_exclusions();
3554
3555        //  Remove temporary structures
3556        delete [] bondsWithAtom;  bondsWithAtom = 0;
3557        delete tmpArena;  tmpArena = 0;
3558
3559        if (exclusions != NULL)
3560       delete [] exclusions;
3561
3562        // 1-4 exclusions which are also fully excluded were eliminated by hash table
3563        numTotalExclusions = exclusionSet.size();
3564        if ( ! CkMyPe() ) {
3565          iout << iINFO << "ADDED " << (numTotalExclusions - numExclusions) << " IMPLICIT EXCLUSIONS\n" << endi;
3566        }
3567        exclusions = new Exclusion[numTotalExclusions];
3568        UniqueSetIter<Exclusion> exclIter(exclusionSet);
3569        for ( exclIter=exclIter.begin(),i=0; exclIter != exclIter.end(); exclIter++,i++ )
3570        {
3571          exclusions[i] = *exclIter;
3572        }
3573        // Free exclusionSet storage
3574        // exclusionSet.clear(1);
3575        exclusionSet.clear();
3576
3577        DebugM(3,"Building exclusion lists.\n");
3578     
3579        for (i=0; i<numAtoms; i++)
3580        {
3581          byAtomSize[i] = 0;
3582        }
3583        numCalcExclusions = 0;
3584        numCalcFullExclusions = 0;
3585        for (i=0; i<numTotalExclusions; i++)
3586        {
3587          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3588                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
3589          byAtomSize[exclusions[i].atom1]++;
3590          numCalcExclusions++;
3591          if ( ! exclusions[i].modified ) numCalcFullExclusions++;
3592        }
3593
3594        for (i=0; i<numAtoms; i++)
3595        {
3596          exclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3597          exclusionsByAtom[i][byAtomSize[i]] = -1;
3598          byAtomSize[i] = 0;
3599        }
3600        for (i=0; i<numTotalExclusions; i++)
3601        {
3602          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3603                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
3604          int a1 = exclusions[i].atom1;
3605          exclusionsByAtom[a1][byAtomSize[a1]++] = i;
3606        }
3607
3608        int32 *byAtomSize2 = new int32[numAtoms];
3609
3610        for (i=0; i<numAtoms; i++)
3611        {
3612          byAtomSize[i] = 0;
3613          byAtomSize2[i] = 0;
3614        }
3615
3616        for (i=0; i<numTotalExclusions; i++)
3617        {
3618          if ( numFixedAtoms && fixedAtomFlags[exclusions[i].atom1]
3619                             && fixedAtomFlags[exclusions[i].atom2] ) continue;
3620          if ( exclusions[i].modified ) {
3621            byAtomSize2[exclusions[i].atom1]++;
3622            byAtomSize2[exclusions[i].atom2]++;
3623          } else {
3624            byAtomSize[exclusions[i].atom1]++;
3625            byAtomSize[exclusions[i].atom2]++;
3626          }
3627        }
3628
3629        for (i=0; i<numAtoms; i++)
3630        {
3631          fullExclusionsByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3632          fullExclusionsByAtom[i][0] = 0;
3633          modExclusionsByAtom[i] = arena->getNewArray(byAtomSize2[i]+1);
3634          modExclusionsByAtom[i][0] = 0;
3635        }
3636
3637        for (i=0; i<numTotalExclusions; i++)
3638        {
3639          int a1 = exclusions[i].atom1;
3640          int a2 = exclusions[i].atom2;
3641          if ( numFixedAtoms && fixedAtomFlags[a1]
3642                             && fixedAtomFlags[a2] ) continue;
3643          int32 *l1, *l2;
3644          if ( exclusions[i].modified ) {
3645            l1 = modExclusionsByAtom[a1];
3646            l2 = modExclusionsByAtom[a2];
3647          } else {
3648            l1 = fullExclusionsByAtom[a1];
3649            l2 = fullExclusionsByAtom[a2];
3650          }
3651          l1[++(*l1)] = a2;
3652          l2[++(*l2)] = a1;
3653        }
3654
3655        if ( ! CkMyPe() && simParams->printExclusions ) {
3656          for (i=0; i<numAtoms; i++) {
3657            int32 *lf = fullExclusionsByAtom[i];
3658            iout << "EXCL " << i << " FULL";
3659            int nf = *(lf++);
3660            for ( int j = 0; j < nf; ++j ) {
3661              iout << " " << *(lf++);
3662            }
3663            iout << "\n";
3664            int32 *lm = modExclusionsByAtom[i];
3665            iout << "EXCL " << i << " MOD";
3666            int nm = *(lm++);
3667            for ( int j = 0; j < nm; ++j ) {
3668              iout << " " << *(lm++);
3669            }
3670            iout << "\n" << endi;
3671          }
3672        }
3673
3674        // DRUDE
3675        if (is_drude_psf || simParams->drudeOn) {
3676
3677          // build Thole (screened Coulomb) correction terms;
3678          // they are constructed implicitly from exclusions
3679
3680          // free the previous Thole array if already allocated
3681          if (tholes != NULL) delete[] tholes;
3682          numTholes = 0;
3683
3684          // count the number of Thole terms
3685          for (i = 0;  i < numTotalExclusions;  i++) {
3686            /* skip over the modified exclusions */
3687            if (exclusions[i].modified) continue;
3688            int a1 = exclusions[i].atom1;
3689            int a2 = exclusions[i].atom2;
3690            if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3691              numTholes++;
3692            }
3693          }
3694
3695          // allocate space for Thole terms
3696          if (numTholes != 0) tholes = new Thole[numTholes];
3697          else tholes = NULL;
3698          int nt = 0;
3699
3700          Real c = COULOMB*simParams->nonbondedScaling/simParams->dielectric;
3701
3702          // store Thole terms
3703          for (i = 0;  i < numTotalExclusions;  i++) {
3704            /* skip over the modified exclusions */
3705            if (exclusions[i].modified) continue;
3706            int a1 = exclusions[i].atom1;
3707            int a2 = exclusions[i].atom2;
3708            // exclusions are stored with a1 < a2
3709            if (a2 < numAtoms-1 && is_drude(a1+1) && is_drude(a2+1)) {
3710              Real thsum = drudeConsts[a1].thole + drudeConsts[a2].thole;
3711              Real aprod = drudeConsts[a1].alpha * drudeConsts[a2].alpha;
3712              // guard against having alpha==0
3713              Real apower = (aprod <= 0 ? 0 : powf(aprod, -1.f/6));
3714              tholes[nt].atom1 = a1;
3715              tholes[nt].atom2 = a1+1;
3716              tholes[nt].atom3 = a2;
3717              tholes[nt].atom4 = a2+1;
3718              tholes[nt].aa = apower * thsum;
3719              tholes[nt].qq = c * atoms[a1+1].charge * atoms[a2+1].charge;
3720              nt++;
3721            }
3722          }
3723
3724          // build Thole lists by atom
3725          DebugM(3, "Building Thole correction term lists.\n");
3726          tholesByAtom = new int32 *[numAtoms];
3727
3728          for (i = 0;  i < numAtoms;  i++) {
3729            byAtomSize[i] = 0;
3730          }
3731          numCalcTholes = 0;
3732          for (i = 0;  i < numTholes;  i++) {
3733            if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3734                               && fixedAtomFlags[tholes[i].atom2]
3735                               && fixedAtomFlags[tholes[i].atom3]
3736                               && fixedAtomFlags[tholes[i].atom4] ) continue;
3737            if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1) continue;
3738            byAtomSize[tholes[i].atom1]++;
3739            numCalcTholes++;
3740          }
3741          for (i = 0;  i < numAtoms;  i++) {
3742            tholesByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3743            tholesByAtom[i][byAtomSize[i]] = -1;
3744            byAtomSize[i] = 0;
3745          }
3746          for (i = 0;  i < numTholes;  i++) {
3747            if ( numFixedAtoms && fixedAtomFlags[tholes[i].atom1]
3748                               && fixedAtomFlags[tholes[i].atom2]
3749                               && fixedAtomFlags[tholes[i].atom3]
3750                               && fixedAtomFlags[tholes[i].atom4] ) continue;
3751            if ( pair_self && fepAtomFlags[tholes[i].atom1] != 1) continue;
3752            int a1 = tholes[i].atom1;
3753            tholesByAtom[a1][byAtomSize[a1]++] = i;
3754          }
3755
3756          // build anisotropic lists by atom
3757          DebugM(3, "Building anisotropic term lists.\n");
3758          anisosByAtom = new int32 *[numAtoms];
3759
3760          for (i = 0;  i < numAtoms;  i++) {
3761            byAtomSize[i] = 0;
3762          }
3763          numCalcAnisos = 0;
3764          for (i = 0;  i < numAnisos;  i++) {
3765            if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
3766                               && fixedAtomFlags[anisos[i].atom2]
3767                               && fixedAtomFlags[anisos[i].atom3]
3768                               && fixedAtomFlags[anisos[i].atom4] ) continue;
3769            if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1) continue;
3770            byAtomSize[anisos[i].atom1]++;
3771            numCalcAnisos++;
3772          }
3773          for (i = 0;  i < numAtoms;  i++) {
3774            anisosByAtom[i] = arena->getNewArray(byAtomSize[i]+1);
3775            anisosByAtom[i][byAtomSize[i]] = -1;
3776            byAtomSize[i] = 0;
3777          }
3778          for (i = 0;  i < numAnisos;  i++) {
3779            if ( numFixedAtoms && fixedAtomFlags[anisos[i].atom1]
3780                               && fixedAtomFlags[anisos[i].atom2]
3781                               && fixedAtomFlags[anisos[i].atom3]
3782                               && fixedAtomFlags[anisos[i].atom4] ) continue;
3783            if ( pair_self && fepAtomFlags[anisos[i].atom1] != 1) continue;
3784            int a1 = anisos[i].atom1;
3785            anisosByAtom[a1][byAtomSize[a1]++] = i;
3786          }
3787
3788        }
3789        // DRUDE
3790
3791        delete [] byAtomSize;  byAtomSize = 0;
3792        delete [] byAtomSize2;  byAtomSize2 = 0;
3793
3794
3795        //  Allocate an array to hold the exclusions for each atom
3796        all_exclusions = new ExclusionCheck[numAtoms];
3797
3798        for (i=0; i<numAtoms; i++)
3799        {
3800          all_exclusions[i].min = numAtoms;
3801          all_exclusions[i].max = -1;
3802        }
3803        for (i=0; i<numTotalExclusions; i++)
3804        {
3805          // first atom should alway have lower number!
3806          int a1 = exclusions[i].atom1;
3807          int a2 = exclusions[i].atom2;
3808          if ( numFixedAtoms && fixedAtomFlags[a1]
3809                             && fixedAtomFlags[a2] ) continue;
3810          if ( all_exclusions[a1].min > a2 ) all_exclusions[a1].min = a2;
3811          if ( all_exclusions[a2].min > a1 ) all_exclusions[a2].min = a1;
3812          if ( a2 > all_exclusions[a1].max ) all_exclusions[a1].max = a2;
3813          if ( a1 > all_exclusions[a2].max ) all_exclusions[a2].max = a1;
3814        }
3815
3816        // build array of all full exclusions for water etc.
3817        int maxDenseAllFull = 0;
3818        int numDenseAllFull = 0;
3819        for (i=0; i<numAtoms; i++) {
3820          int iInMiddle = ( i < all_exclusions[i].max &&
3821                        i > all_exclusions[i].min ) ? 1 : 0;
3822          int s = all_exclusions[i].max - all_exclusions[i].min + 1;
3823          if ( s == fullExclusionsByAtom[i][0] + iInMiddle ) {
3824             if ( s > maxDenseAllFull ) maxDenseAllFull = s;
3825             all_exclusions[i].flags = (char*)-1;  // shared array
3826          } else {
3827             all_exclusions[i].flags = 0;  // individual array
3828          }
3829        }
3830        char *denseFullArray = exclArena->getNewArray(maxDenseAllFull);
3831        for ( i=0; i<maxDenseAllFull; ++i ) denseFullArray[i] = EXCHCK_FULL;
3832
3833        int exclmem = maxDenseAllFull;
3834        int maxExclusionFlags = simParams->maxExclusionFlags;
3835        for (i=0; i<numAtoms; i++) {
3836          int s = all_exclusions[i].max - all_exclusions[i].min + 1;
3837          if ( all_exclusions[i].max != -1 ) {
3838            if ( all_exclusions[i].flags ) {
3839              all_exclusions[i].flags = denseFullArray;
3840              ++numDenseAllFull;
3841            } else if ( s < maxExclusionFlags ) {
3842              char *f = all_exclusions[i].flags = exclArena->getNewArray(s);
3843              for ( int k=0; k<s; ++k ) f[k] = 0;
3844              exclmem += s;
3845            } else {
3846              all_exclusions[i].flags = 0;  // need to build on the fly
3847            }
3848          } else {
3849            all_exclusions[i].flags = (char*)-1; // should never dereference
3850          }
3851        }
3852        if ( 0 ) {
3853          iout << iINFO << numTotalExclusions << " exclusions consume "
3854             << exclmem << " bytes.\n" << endi;
3855          iout << iINFO << numDenseAllFull
3856             << " atoms sharing one array.\n" << endi;
3857        }
3858        for (i=0; i<numTotalExclusions; i++)
3859        {
3860          int a1 = exclusions[i].atom1;
3861          int a2 = exclusions[i].atom2;
3862          if ( numFixedAtoms && fixedAtomFlags[a1]
3863                             && fixedAtomFlags[a2] ) continue;
3864          if ( exclusions[i].modified ) {
3865            if ( all_exclusions[a1].flags )
3866              all_exclusions[a1].flags[a2-all_exclusions[a1].min] = EXCHCK_MOD;
3867            if ( all_exclusions[a2].flags )
3868              all_exclusions[a2].flags[a1-all_exclusions[a2].min] = EXCHCK_MOD;
3869          } else {
3870            if ( all_exclusions[a1].flags )
3871              all_exclusions[a1].flags[a2-all_exclusions[a1].min] = EXCHCK_FULL;
3872            if ( all_exclusions[a2].flags )
3873              all_exclusions[a2].flags[a1-all_exclusions[a2].min] = EXCHCK_FULL;
3874          }
3875        }
3876     }
3877
3878     /*    END OF FUNCTION build_lists_by_atom    */
3879
3880     /****************************************************************/
3881     /*                */
3882     /*      FUNCTION build_exclusions    */
3883     /*                */
3884     /*  This function builds a list of all the exlcusions       */
3885     /*  atoms.  These lists include explicit exclusions as well as  */
3886     /*  exclusions that are calculated based on the bonded structure*/
3887     /*  and the exclusion flag.  For each pair of atoms that are    */
3888     /*  excluded, the larger of the 2 atom indexes is stored in the */
3889     /*  array of the smaller index.  All the arrays are not sorted. */
3890     /*  Then to determine if two atoms have an exclusion, a linear  */
3891     /*  search is done on the array of the atom with the smaller    */
3892     /*  index for the larger index.          */
3893     /*  If the exclusion policy is set to scaled1-4, there are  */
3894     /*  actually two lists built.  One contains the pairs of atoms  */
3895     /*  that are to be exlcuded (i.e., explicit exclusions, 1-2,    */
3896     /*  and 1-3 interactions) and the other contains just the 1-4   */
3897     /*  interactions, since they will need to be determined   */
3898     /*  independantly of the other exclusions.      */
3899     /*                */
3900     /****************************************************************/
3901
3902     void Molecule::build_exclusions()
3903     {
3904       register int i;          //  Loop counter
3905       ExclusionSettings exclude_flag;    //  Exclusion policy
3906
3907       exclude_flag = simParams->exclude;
3908
3909       //  Go through the explicit exclusions and add them to the arrays
3910       for (i=0; i<numExclusions; i++)
3911       {
3912         exclusionSet.add(exclusions[i]);
3913       }
3914
3915       // If this is AMBER force field, and readExclusions is TRUE,
3916       // then all the exclusions were read from parm file, and we
3917       // shouldn't generate any of them.
3918       if (!simParams->amberOn || !simParams->readExclusions)
3919       { //  Now calculate the bonded exlcusions based on the exclusion policy
3920         switch (exclude_flag)
3921         {
3922          case NONE:
3923            break;
3924          case ONETWO:
3925            build12excl();
3926            break;
3927           case ONETHREE:
3928             build12excl();
3929             build13excl();
3930             break;
3931           case ONEFOUR:
3932             build12excl();
3933             build13excl();
3934             build14excl(0);
3935             break;
3936           case SCALED14:
3937             build12excl();
3938             build13excl();
3939             build14excl(1);
3940             break;
3941         }
3942       }
3943
3944       stripFepExcl();
3945
3946       // DRUDE
3947       if (is_lonepairs_psf || is_drude_psf) {
3948         build_inherited_excl(SCALED14 == exclude_flag);
3949       }
3950     }
3951     /*      END OF FUNCTION build_exclusions    */
3952
3953
3954     // Extend exclusions for the Drude model.  The Drude model is generally
3955     // used with the 1-3 exclusion policy, although the code below also
3956     // supports the 1-2 exclusion policy.  The use of light (or massless)
3957     // pseudo-atoms requires the introduction of extra exclusions.
3958     //
3959     // Here is the algorithm for determining Drude model exclusions:
3960     // (1)  Each Drude particle and each lone pair has a single parent atom.
3961     //      The parent atom must be a heavy atom.
3962     // (2)  Each Drude particle and lone pair inherit the exclusions of its
3963     //      parent atom.
3964     // (3)  If two heavy atoms are excluded and they both have either a
3965     //      Drude particle or a lone pair, the these light (or massless)
3966     //      particles are also excluded from interacting with each other.
3967     void Molecule::build_inherited_excl(int modified) {
3968       ExclusionSettings exclude_flag = simParams->exclude;
3969       int32 *bond1, *bond2, *bond3, *bond4, *bond5;
3970       int32 i, j, mid1, mid2, mid3, mid4;
3971
3972       // validate that each Drude or lone pair particle
3973       // has a unique parent that is a heavy atom
3974       for (i = 0;  i < numAtoms;  i++) {
3975
3976         if (!is_drude(i) && !is_lp(i)) continue;
3977         // make sure that i is either Drude or LP
3978
3979         // find parent (heavy) atom of particle i
3980         bond1 = bondsWithAtom[i];
3981
3982         if (-1 == *bond1) {  // i must have one bond
3983           char err_msg[512];
3984           const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
3985           sprintf(err_msg, "FOUND ISOLATED %s PARTICLE %d", idescrip, i+1);
3986           NAMD_die(err_msg);
3987         }
3988         if (-1 != *(bond1+1)) {  // and only one bond
3989           char err_msg[512];
3990           const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
3991           sprintf(err_msg, "FOUND MULTIPLY LINKED %s PARTICLE %d",
3992               idescrip, i+1);
3993           NAMD_die(err_msg);
3994         }
3995
3996         // mid1 is parent of particle i
3997         mid1 = bonds[*bond1].atom1;
3998         if (mid1 == i) mid1 = bonds[*bond1].atom2;
3999
4000         // make sure that mid1 is a heavy atom
4001         if (is_drude(mid1) || is_lp(mid1) || is_hydrogen(mid1)) {
4002           char err_msg[512];
4003           const char *idescrip = (is_drude(i) ? "DRUDE" : "LONE PAIR");
4004           sprintf(err_msg, "PARENT ATOM %d of %s PARTICLE %d "
4005               "IS NOT HEAVY ATOM", mid1+1, idescrip, i+1);
4006           NAMD_die(err_msg);
4007         }
4008
4009         if (exclude_flag == NONE) {
4010           // add (i,mid1) as an exclusion
4011           if (i < mid1) {
4012             exclusionSet.add(Exclusion(i, mid1));
4013           }
4014           else {
4015             exclusionSet.add(Exclusion(mid1, i));
4016           }
4017
4018           // also exclude any Drude particles or LPs bonded to mid1
4019           bond2 = bondsWithAtom[mid1];
4020           while (*bond2 != -1) {
4021             j = bonds[*bond2].atom1;
4022             if ((is_drude(j) || is_lp(j)) && j != mid1) {
4023               if      (i < j) exclusionSet.add(Exclusion(i, j));
4024               else if (j < i) exclusionSet.add(Exclusion(j, i));
4025             }
4026             j = bonds[*bond2].atom2;
4027             if ((is_drude(j) || is_lp(j)) && j != mid1) {
4028               if      (i < j) exclusionSet.add(Exclusion(i, j));
4029               else if (j < i) exclusionSet.add(Exclusion(j, i));
4030             }
4031             bond2++;
4032           }
4033         }
4034         else {  // if ONETWO or ONETHREE or ONEFOUR or SCALED14
4035
4036           // find the next link
4037           bond2 = bondsWithAtom[mid1];
4038
4039           // loop through all the bonds connected to atom mid1
4040           while (*bond2 != -1) {
4041             if (bonds[*bond2].atom1 == mid1) {
4042               mid2 = bonds[*bond2].atom2;
4043             }
4044             else {
4045               mid2 = bonds[*bond2].atom1;
4046             }
4047
4048             // Make sure that we don't double back to where we started from.
4049             // Doing so causes strange behavior.
4050             if (mid2 == i) {
4051               bond2++;
4052               continue;
4053             }
4054
4055             if (exclude_flag == ONETWO) {
4056               // add (i,mid2) as an exclusion
4057               if (i < mid2) {
4058                 exclusionSet.add(Exclusion(i, mid2));
4059               }
4060               else {
4061                 exclusionSet.add(Exclusion(mid2, i));
4062               }
4063
4064               // also exclude any Drude particles or LPs bonded to mid2
4065               bond3 = bondsWithAtom[mid2];
4066               while (*bond3 != -1) {
4067                 j = bonds[*bond3].atom1;
4068                 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4069                   if      (i < j) exclusionSet.add(Exclusion(i, j));
4070                   else if (j < i) exclusionSet.add(Exclusion(j, i));
4071                 }
4072                 j = bonds[*bond3].atom2;
4073                 if ((is_drude(j) || is_lp(j)) && j != mid2) {
4074                   if      (i < j) exclusionSet.add(Exclusion(i, j));
4075                   else if (j < i) exclusionSet.add(Exclusion(j, i));
4076                 }
4077                 bond3++;
4078               }
4079             }
4080             else { // if ONETHREE or ONEFOUR or SCALED14
4081
4082               // find the next link
4083               bond3 = bondsWithAtom[mid2];
4084
4085               // loop through all the bonds connected to mid2
4086               while (*bond3 != -1) {
4087
4088                 if (bonds[*bond3].atom1 == mid2) {
4089                   mid3 = bonds[*bond3].atom2;
4090                 }
4091                 else {
4092                   mid3 = bonds[*bond3].atom1;
4093                 }
4094
4095                 // Make sure we don't double back to where we started.
4096                 // Doing so causes strange behavior.
4097                 if (mid3 == mid1) {
4098                   bond3++;
4099                   continue;
4100                 }
4101
4102                 // add (i,mid3) as an exclusion
4103                 if (i < mid3) {
4104                   exclusionSet.add(Exclusion(i, mid3));
4105                 }
4106                 else if (mid3 < i) {
4107                   exclusionSet.add(Exclusion(mid3, i));
4108                 }
4109
4110                 if (exclude_flag == ONETHREE) {
4111                   // also exclude any Drude particles or LPs bonded to mid3
4112                   bond4 = bondsWithAtom[mid3];
4113                   while (*bond4 != -1) {
4114                     j = bonds[*bond4].atom1;
4115                     if ((is_drude(j) || is_lp(j)) && j != mid3) {
4116                       if      (i < j) exclusionSet.add(Exclusion(i, j));
4117                       else if (j < i) exclusionSet.add(Exclusion(j, i));
4118                     }
4119                     j = bonds[*bond4].atom2;
4120                     if ((is_drude(j) || is_lp(j)) && j != mid3) {
4121                       if      (i < j) exclusionSet.add(Exclusion(i, j));
4122                       else if (j < i) exclusionSet.add(Exclusion(j, i));
4123                     }
4124                     bond4++;
4125                   }
4126                 }
4127                 else { // if ONEFOUR or SCALED14
4128
4129                   // find next link
4130                   bond4 = bondsWithAtom[mid3];
4131
4132                   // loop through all the bonds connected to mid3
4133                   while (*bond4 != -1) {
4134
4135                     if (bonds[*bond4].atom1 == mid3) {
4136                       mid4 = bonds[*bond4].atom2;
4137                     }
4138                     else {
4139                       mid4 = bonds[*bond4].atom1;
4140                     }
4141
4142                     // Make sure we don't double back to where we started.
4143                     // Doing so causes strange behavior.
4144                     if (mid4 == mid2) {
4145                       bond4++;
4146                       continue;
4147                     }
4148
4149                     if (is_drude(mid4) || is_lp(mid4)) {
4150                       // (i,mid4) is 1-3 excl
4151                       if (i < mid4) {
4152                         exclusionSet.add(Exclusion(i, mid4));
4153                       }
4154                       else if (mid4 < i) {
4155                         exclusionSet.add(Exclusion(mid4, i));
4156                       }
4157                       bond4++;
4158                       continue;
4159                     }
4160
4161                     // (mid1,mid4) is an existing heavy atom exclusion
4162                     // if we have modified 1-4 exclusions, make sure
4163                     // that (mid1,mid4) is modified 1-4 exclusion
4164                     // rather than something closer due to a ring
4165                     int modi = modified;
4166                     if (modified) {
4167                       int amin = (mid1 < mid4 ? mid1 : mid4);
4168                       int amax = (mid1 >= mid4 ? mid1 : mid4);
4169                       Exclusion *pe = exclusionSet.find(Exclusion(amin,amax));
4170                       if (pe==0) {
4171                         // since there is not an existing exclusion
4172                         // between (mid1,mid4), don't inherit!
4173                         bond4++;
4174                         continue;
4175                       }
4176                       modi = pe->modified;
4177                     }
4178
4179                     if (i < mid4) {
4180                       exclusionSet.add(Exclusion(i, mid4, modi));
4181                     }
4182                     else if (mid4 < i) {
4183                       exclusionSet.add(Exclusion(mid4, i, modi));                    
4184                     }
4185
4186                     // also exclude any Drude particles or LPs bonded to mid4
4187                     // using the "modi" setting of (mid1,mid4) exclusion
4188                     bond5 = bondsWithAtom[mid4];
4189                     while (*bond5 != -1) {
4190                       j = bonds[*bond5].atom1;
4191                       if ((is_drude(j) || is_lp(j)) && j != mid4) {
4192                         if      (i<j) exclusionSet.add(Exclusion(i,j,modi));
4193                         else if (j<i) exclusionSet.add(Exclusion(j,i,modi));
4194                       }
4195                       j = bonds[*bond5].atom2;
4196                       if ((is_drude(j) || is_lp(j)) && j != mid4) {
4197                         if      (i<j) exclusionSet.add(Exclusion(i,j,modi));
4198                         else if (j<i) exclusionSet.add(Exclusion(j,i,modi));
4199                       }
4200                       bond5++;
4201                     }
4202                     ++bond4;
4203                   } // while bond4
4204
4205                 } // else (if ONEFOUR or SCALED14)
4206
4207                 ++bond3;
4208               } // while bond3
4209
4210             } // else (if ONETHREE or ONEFOUR or SCALED14)
4211            
4212             ++bond2;
4213           } // while bond2
4214
4215         } // else (if ONETWO or ONETHREE or ONEFOUR or SCALED14)
4216
4217       } // for i
4218     } 
4219     // DRUDE
4220
4221
4222     /************************************************************************/
4223     /*                  */
4224     /*      FUNCTION build12excl        */
4225     /*                  */
4226     /************************************************************************/
4227
4228     void Molecule::build12excl(void)       
4229     {
4230        int32 *current_val;  //  Current value to check
4231        register int i;    //  Loop counter to loop through all atoms
4232        
4233        //  Loop through all the atoms marking the bonded interactions for each one
4234        for (i=0; i<numAtoms; i++)
4235        {
4236       current_val = bondsWithAtom[i];
4237        
4238       //  Loop through all the bonds for this atom
4239       while (*current_val != -1)
4240       {
4241          if (bonds[*current_val].atom1 == i)
4242          {
4243       if (i<bonds[*current_val].atom2)
4244       {
4245          exclusionSet.add(Exclusion(i,bonds[*current_val].atom2));
4246       }
4247          }
4248          else
4249          {
4250       if (i<bonds[*current_val].atom1)
4251       {
4252          exclusionSet.add(Exclusion(i,bonds[*current_val].atom1));
4253       }
4254          }
4255     
4256          ++current_val;