OOPSE/libmdtools/ForceFields.cpp

#include <cstdlib>

#ifdef IS_MPI
#include <mpi.h>
#include <mpi++.h>
#endif // is_mpi


#include "simError.h"
#include "ForceFields.hpp"
#include "Atom.hpp"
#include "fortranWrappers.hpp"


void ForceFields::calcRcut( void ){

#ifdef IS_MPI
  double tempBig = bigSigma;
  MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
#endif  //is_mpi

  //calc rCut and rList

  entry_plug->rCut = 2.5 * bigSigma;
  if(entry_plug->rCut > (entry_plug->box_x / 2.0)) 
    entry_plug->rCut = entry_plug->box_x / 2.0;
  if(entry_plug->rCut > (entry_plug->box_y / 2.0)) 
    entry_plug->rCut = entry_plug->box_y / 2.0;
  if(entry_plug->rCut > (entry_plug->box_z / 2.0)) 
    entry_plug->rCut = entry_plug->box_z / 2.0;
  
  entry_plug->rList = entry_plug->rCut + 1.0;
  
}

void ForceFields::doForces( int calcPot, int calcStress ){

  int i, isError;
  double* frc;
  double* pos;
  double* trq;
  double* tau;
  double* A;
  double* u_l;;
  DirectionalAtom* dAtom;

  double ut[3];

  //u_l = new double[entry_plug->n_atoms];

  short int passedCalcPot = (short int)calcPot;
  short int passedCalcStress = (short int)calcStress;

  // forces are zeroed here, before any are acumulated.
  // NOTE: do not rezero the forces in Fortran.


  for(i=0; i<entry_plug->n_atoms; i++){
    entry_plug->atoms[i]->zeroForces();

    if( entry_plug->atoms[i]->isDirectional() ){
      dAtom = (DirectionalAtom *)entry_plug->atoms[i];
      dAtom->getU(ut);
    

      if(dAtom->getIndex()== 1){
        std::cerr << "atom 2's u_l = " << ut[0] << ", " << ut[1]
                  << ", " << ut[2] << "\n";
      }
    }
    
  }

  for(i=0; i<entry_plug->n_mol; i++ ){
    entry_plug->molecules[i].calcForces();
  }

  frc = Atom::getFrcArray();
  pos = Atom::getPosArray();
  trq = Atom::getTrqArray();
  A   = Atom::getAmatArray();
  u_l = Atom::getUlArray();
  tau = entry_plug->tau;

  
  isError = 0;
  entry_plug->lrPot = 0.0;


  fortranForceLoop( pos,
                    A,
                    u_l,
                    frc,
                    trq,
                    tau,
                    &(entry_plug->lrPot), 
                    &passedCalcPot,
                    &passedCalcStress,
                    &isError );


  //  delete[] u_l;

  if( isError ){
    sprintf( painCave.errMsg,
             "Error returned from the fortran force calculation.\n" );
    painCave.isFatal = 1;
    simError();
  }

#ifdef IS_MPI
  sprintf( checkPointMsg,
           "returned from the force calculation.\n" );
  MPIcheckPoint();
#endif // is_mpi

}


void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
  
  int isError;
  
  isError = 0;
  initFortranFF( &ljMixPolicy, &useReactionField, &isError );
  
  if(isError){
    sprintf( painCave.errMsg,
             "ForceField error: There was an error initializing the forceField in fortran.\n" );
    painCave.isFatal = 1;
    simError();
  }

  
#ifdef IS_MPI
  sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
  MPIcheckPoint();
#endif // is_mpi
  
}
Revision:	443
Committed:	Wed Apr 2 22:19:03 2003 UTC (22 years, 7 months ago) by mmeineke
File size:	2904 byte(s)
Log Message:	dipoles mostly work, but there is a memory leak somewhere.
#	User	Rev	Content
1	mmeineke	377	#include <cstdlib>
2
3			#ifdef IS_MPI
4			#include <mpi.h>
5	mmeineke	428	#include <mpi++.h>
6	mmeineke	377	#endif // is_mpi
7
8
9			#include "simError.h"
10			#include "ForceFields.hpp"
11			#include "Atom.hpp"
12			#include "fortranWrappers.hpp"
13
14
15	mmeineke	420	void ForceFields::calcRcut( void ){
16
17			#ifdef IS_MPI
18			double tempBig = bigSigma;
19			MPI::COMM_WORLD.Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX );
20			#endif //is_mpi
21
22			//calc rCut and rList
23
24	mmeineke	428	entry_plug->rCut = 2.5 * bigSigma;
25	mmeineke	420	if(entry_plug->rCut > (entry_plug->box_x / 2.0))
26			entry_plug->rCut = entry_plug->box_x / 2.0;
27			if(entry_plug->rCut > (entry_plug->box_y / 2.0))
28			entry_plug->rCut = entry_plug->box_y / 2.0;
29			if(entry_plug->rCut > (entry_plug->box_z / 2.0))
30			entry_plug->rCut = entry_plug->box_z / 2.0;
31
32			entry_plug->rList = entry_plug->rCut + 1.0;
33
34			}
35
36	mmeineke	377	void ForceFields::doForces( int calcPot, int calcStress ){
37
38			int i, isError;
39			double* frc;
40			double* pos;
41			double* trq;
42			double* tau;
43			double* A;
44	mmeineke	443	double* u_l;;
45			DirectionalAtom* dAtom;
46	mmeineke	377
47	mmeineke	443	double ut[3];
48
49			//u_l = new double[entry_plug->n_atoms];
50
51	mmeineke	377	short int passedCalcPot = (short int)calcPot;
52			short int passedCalcStress = (short int)calcStress;
53
54			// forces are zeroed here, before any are acumulated.
55			// NOTE: do not rezero the forces in Fortran.
56
57	mmeineke	443
58	mmeineke	377	for(i=0; i<entry_plug->n_atoms; i++){
59			entry_plug->atoms[i]->zeroForces();
60	mmeineke	443
61			if( entry_plug->atoms[i]->isDirectional() ){
62			dAtom = (DirectionalAtom *)entry_plug->atoms[i];
63			dAtom->getU(ut);
64
65
66			if(dAtom->getIndex()== 1){
67			std::cerr << "atom 2's u_l = " << ut[0] << ", " << ut[1]
68			<< ", " << ut[2] << "\n";
69			}
70			}
71
72	mmeineke	377	}
73
74	mmeineke	423	for(i=0; i<entry_plug->n_mol; i++ ){
75	mmeineke	428	entry_plug->molecules[i].calcForces();
76	mmeineke	389	}
77
78	mmeineke	377	frc = Atom::getFrcArray();
79			pos = Atom::getPosArray();
80			trq = Atom::getTrqArray();
81			A = Atom::getAmatArray();
82			u_l = Atom::getUlArray();
83	mmeineke	424	tau = entry_plug->tau;
84	mmeineke	377
85	mmeineke	424
86	mmeineke	377	isError = 0;
87			entry_plug->lrPot = 0.0;
88	mmeineke	393
89	chuckv	438
90	mmeineke	393
91	mmeineke	377	fortranForceLoop( pos,
92			A,
93			u_l,
94			frc,
95			trq,
96			tau,
97			&(entry_plug->lrPot),
98			&passedCalcPot,
99			&passedCalcStress,
100			&isError );
101
102
103	mmeineke	443	// delete[] u_l;
104	chuckv	438
105	mmeineke	377	if( isError ){
106			sprintf( painCave.errMsg,
107			"Error returned from the fortran force calculation.\n" );
108			painCave.isFatal = 1;
109			simError();
110			}
111
112			#ifdef IS_MPI
113			sprintf( checkPointMsg,
114			"returned from the force calculation.\n" );
115			MPIcheckPoint();
116			#endif // is_mpi
117
118			}
119
120
121			void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
122
123			int isError;
124
125			isError = 0;
126			initFortranFF( &ljMixPolicy, &useReactionField, &isError );
127
128			if(isError){
129			sprintf( painCave.errMsg,
130			"ForceField error: There was an error initializing the forceField in fortran.\n" );
131			painCave.isFatal = 1;
132			simError();
133			}
134
135
136			#ifdef IS_MPI
137			sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
138			MPIcheckPoint();
139			#endif // is_mpi
140
141			}