ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/ForceFields.cpp
Revision: 597
Committed: Mon Jul 14 21:28:54 2003 UTC (20 years, 11 months ago) by mmeineke
File size: 2615 byte(s)
Log Message:
found a bug. Unit vectors were not being updated

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 597
66    
67 mmeineke 377 frc = Atom::getFrcArray();
68     pos = Atom::getPosArray();
69     trq = Atom::getTrqArray();
70     A = Atom::getAmatArray();
71     u_l = Atom::getUlArray();
72 mmeineke 561
73 mmeineke 377 isError = 0;
74     entry_plug->lrPot = 0.0;
75 mmeineke 393
76 chuckv 479 for (i=0; i<9; i++) {
77     entry_plug->tau[i] = 0.0;
78     }
79 chuckv 438
80 mmeineke 377 fortranForceLoop( pos,
81     A,
82     u_l,
83     frc,
84     trq,
85 chuckv 479 entry_plug->tau,
86 mmeineke 377 &(entry_plug->lrPot),
87     &passedCalcPot,
88     &passedCalcStress,
89     &isError );
90    
91     if( isError ){
92     sprintf( painCave.errMsg,
93     "Error returned from the fortran force calculation.\n" );
94     painCave.isFatal = 1;
95     simError();
96     }
97    
98     #ifdef IS_MPI
99     sprintf( checkPointMsg,
100     "returned from the force calculation.\n" );
101     MPIcheckPoint();
102     #endif // is_mpi
103 mmeineke 597
104 mmeineke 377
105     }
106    
107    
108     void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
109    
110     int isError;
111    
112     isError = 0;
113     initFortranFF( &ljMixPolicy, &useReactionField, &isError );
114    
115     if(isError){
116     sprintf( painCave.errMsg,
117     "ForceField error: There was an error initializing the forceField in fortran.\n" );
118     painCave.isFatal = 1;
119     simError();
120     }
121    
122    
123     #ifdef IS_MPI
124     sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
125     MPIcheckPoint();
126     #endif // is_mpi
127    
128     }