ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/ForceFields.cpp
Revision: 572
Committed: Wed Jul 2 21:26:55 2003 UTC (21 years ago) by mmeineke
File size: 2606 byte(s)
Log Message:
fixed the bugs introduced by switching the periodic box to a matrix

File Contents

# User Rev Content
1 mmeineke 561 #include <iostream>
2    
3     using namespace std;
4    
5    
6 mmeineke 377 #include <cstdlib>
7    
8     #ifdef IS_MPI
9     #include <mpi.h>
10     #endif // is_mpi
11    
12    
13     #include "simError.h"
14     #include "ForceFields.hpp"
15     #include "Atom.hpp"
16     #include "fortranWrappers.hpp"
17    
18    
19 mmeineke 420 void ForceFields::calcRcut( void ){
20    
21     #ifdef IS_MPI
22     double tempBig = bigSigma;
23 mmeineke 447 MPI_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX,
24     MPI_COMM_WORLD);
25 mmeineke 420 #endif //is_mpi
26    
27     //calc rCut and rList
28    
29 mmeineke 428 entry_plug->rCut = 2.5 * bigSigma;
30 mmeineke 572 if(entry_plug->rCut > (entry_plug->boxLx / 2.0))
31     entry_plug->rCut = entry_plug->boxLx / 2.0;
32     if(entry_plug->rCut > (entry_plug->boxLy / 2.0))
33     entry_plug->rCut = entry_plug->boxLy / 2.0;
34     if(entry_plug->rCut > (entry_plug->boxLz / 2.0))
35     entry_plug->rCut = entry_plug->boxLz / 2.0;
36 mmeineke 420
37     entry_plug->rList = entry_plug->rCut + 1.0;
38    
39     }
40    
41 mmeineke 377 void ForceFields::doForces( int calcPot, int calcStress ){
42    
43     int i, isError;
44     double* frc;
45     double* pos;
46     double* trq;
47     double* A;
48 mmeineke 443 double* u_l;;
49     DirectionalAtom* dAtom;
50 mmeineke 377
51     short int passedCalcPot = (short int)calcPot;
52     short int passedCalcStress = (short int)calcStress;
53    
54 gezelter 484 // forces are zeroed here, before any are accumulated.
55 mmeineke 377 // NOTE: do not rezero the forces in Fortran.
56    
57     for(i=0; i<entry_plug->n_atoms; i++){
58 gezelter 484 entry_plug->atoms[i]->zeroForces();
59 mmeineke 377 }
60    
61 mmeineke 423 for(i=0; i<entry_plug->n_mol; i++ ){
62 mmeineke 428 entry_plug->molecules[i].calcForces();
63 mmeineke 389 }
64    
65 mmeineke 377 frc = Atom::getFrcArray();
66     pos = Atom::getPosArray();
67     trq = Atom::getTrqArray();
68     A = Atom::getAmatArray();
69     u_l = Atom::getUlArray();
70 mmeineke 561
71 mmeineke 377 isError = 0;
72     entry_plug->lrPot = 0.0;
73 mmeineke 393
74 chuckv 479 for (i=0; i<9; i++) {
75     entry_plug->tau[i] = 0.0;
76     }
77 chuckv 438
78 mmeineke 377 fortranForceLoop( pos,
79     A,
80     u_l,
81     frc,
82     trq,
83 chuckv 479 entry_plug->tau,
84 mmeineke 377 &(entry_plug->lrPot),
85     &passedCalcPot,
86     &passedCalcStress,
87     &isError );
88    
89     if( isError ){
90     sprintf( painCave.errMsg,
91     "Error returned from the fortran force calculation.\n" );
92     painCave.isFatal = 1;
93     simError();
94     }
95    
96     #ifdef IS_MPI
97     sprintf( checkPointMsg,
98     "returned from the force calculation.\n" );
99     MPIcheckPoint();
100     #endif // is_mpi
101    
102     }
103    
104    
105     void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
106    
107     int isError;
108    
109     isError = 0;
110     initFortranFF( &ljMixPolicy, &useReactionField, &isError );
111    
112     if(isError){
113     sprintf( painCave.errMsg,
114     "ForceField error: There was an error initializing the forceField in fortran.\n" );
115     painCave.isFatal = 1;
116     simError();
117     }
118    
119    
120     #ifdef IS_MPI
121     sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
122     MPIcheckPoint();
123     #endif // is_mpi
124    
125     }