OOPSE/libmdtools/ForceFields.cpp

#include <iostream>

using namespace std;


#include <stdlib.h>

#ifdef IS_MPI
#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_Allreduce( &tempBig, &bigSigma, 1, MPI_DOUBLE, MPI_MAX,
                 MPI_COMM_WORLD);
#endif  //is_mpi

  //calc rCut and rList

  entry_plug->setDefaultRcut( 2.5 * bigSigma );  
    
}

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

  int i, isError;
  double* frc;
  double* pos;
  double* trq;
  double* A;
  double* u_l;;
  SimState* config;

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

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

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

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

  config = entry_plug->getConfiguration();
  
  frc = config->getFrcArray();
  pos = config->getPosArray();
  trq = config->getTrqArray();
  A   = config->getAmatArray();
  u_l = config->getUlArray();

  isError = 0;
  entry_plug->lrPot = 0.0;

  for (i=0; i<9; i++) {
    entry_plug->tau[i] = 0.0;
  }

  std::cerr << "before\n"
            << "  force[0] = " << frc[0] << "\n"
            << "  pos[0]   = " << pos[0] << "\n"
            << "  trq[0]   = " << trq[0] << "\n"
            << "  A[0]     = " << A[0] << "\n"
            << "  tau      = " << entry_plug->tau[0] << ", " << entry_plug->tau[1] << ", " << entry_plug->tau[2] << "\n"
            << "             " << entry_plug->tau[3] << ", " << entry_plug->tau[4] << ", " << entry_plug->tau[5] << "\n"
            << "             " << entry_plug->tau[6] << ", " << entry_plug->tau[7] << ", " << entry_plug->tau[8] << "\n\n"; 

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


  std::cerr << "after\n"
            << "  force[0] = " << frc[0] << "\n"
            << "  pos[0]   = " << pos[0] << "\n"
            << "  trq[0]   = " << trq[0] << "\n"
            << "  A[0]     = " << A[0] << "\n"
            << "  tau      = " << entry_plug->tau[0] << ", " << entry_plug->tau[1] << ", " << entry_plug->tau[2] << "\n"
            << "             " << entry_plug->tau[3] << ", " << entry_plug->tau[4] << ", " << entry_plug->tau[5] << "\n"
            << "             " << entry_plug->tau[6] << ", " << entry_plug->tau[7] << ", " << entry_plug->tau[8] << "\n\n";

  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:	841
Committed:	Wed Oct 29 17:55:28 2003 UTC (20 years, 8 months ago) by mmeineke
File size:	3403 byte(s)
Log Message:	som efixes to the way rcut is setup, as well as additional debugging comments.
#	User	Rev	Content
1	mmeineke	561	#include <iostream>
2
3			using namespace std;
4
5
6	gezelter	829	#include <stdlib.h>
7	mmeineke	377
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	841	entry_plug->setDefaultRcut( 2.5 * bigSigma );
30	tim	781
31	mmeineke	420	}
32
33	mmeineke	377	void ForceFields::doForces( int calcPot, int calcStress ){
34
35			int i, isError;
36			double* frc;
37			double* pos;
38			double* trq;
39			double* A;
40	mmeineke	443	double* u_l;;
41	mmeineke	670	SimState* config;
42	mmeineke	377
43			short int passedCalcPot = (short int)calcPot;
44			short int passedCalcStress = (short int)calcStress;
45
46	gezelter	484	// forces are zeroed here, before any are accumulated.
47	mmeineke	377	// NOTE: do not rezero the forces in Fortran.
48
49			for(i=0; i<entry_plug->n_atoms; i++){
50	gezelter	484	entry_plug->atoms[i]->zeroForces();
51	mmeineke	377	}
52
53	mmeineke	423	for(i=0; i<entry_plug->n_mol; i++ ){
54	mmeineke	428	entry_plug->molecules[i].calcForces();
55	mmeineke	389	}
56
57	mmeineke	670	config = entry_plug->getConfiguration();
58	mmeineke	597
59	mmeineke	670	frc = config->getFrcArray();
60			pos = config->getPosArray();
61			trq = config->getTrqArray();
62			A = config->getAmatArray();
63			u_l = config->getUlArray();
64	mmeineke	561
65	mmeineke	377	isError = 0;
66			entry_plug->lrPot = 0.0;
67	mmeineke	393
68	chuckv	479	for (i=0; i<9; i++) {
69			entry_plug->tau[i] = 0.0;
70			}
71	chuckv	438
72	mmeineke	841	std::cerr << "before\n"
73			<< " force[0] = " << frc[0] << "\n"
74			<< " pos[0] = " << pos[0] << "\n"
75			<< " trq[0] = " << trq[0] << "\n"
76			<< " A[0] = " << A[0] << "\n"
77			<< " tau = " << entry_plug->tau[0] << ", " << entry_plug->tau[1] << ", " << entry_plug->tau[2] << "\n"
78			<< " " << entry_plug->tau[3] << ", " << entry_plug->tau[4] << ", " << entry_plug->tau[5] << "\n"
79			<< " " << entry_plug->tau[6] << ", " << entry_plug->tau[7] << ", " << entry_plug->tau[8] << "\n\n";
80
81	mmeineke	377	fortranForceLoop( pos,
82			A,
83			u_l,
84			frc,
85			trq,
86	chuckv	479	entry_plug->tau,
87	mmeineke	377	&(entry_plug->lrPot),
88			&passedCalcPot,
89			&passedCalcStress,
90			&isError );
91
92	mmeineke	841
93			std::cerr << "after\n"
94			<< " force[0] = " << frc[0] << "\n"
95			<< " pos[0] = " << pos[0] << "\n"
96			<< " trq[0] = " << trq[0] << "\n"
97			<< " A[0] = " << A[0] << "\n"
98			<< " tau = " << entry_plug->tau[0] << ", " << entry_plug->tau[1] << ", " << entry_plug->tau[2] << "\n"
99			<< " " << entry_plug->tau[3] << ", " << entry_plug->tau[4] << ", " << entry_plug->tau[5] << "\n"
100			<< " " << entry_plug->tau[6] << ", " << entry_plug->tau[7] << ", " << entry_plug->tau[8] << "\n\n";
101
102	mmeineke	377	if( isError ){
103			sprintf( painCave.errMsg,
104			"Error returned from the fortran force calculation.\n" );
105			painCave.isFatal = 1;
106			simError();
107			}
108
109			#ifdef IS_MPI
110			sprintf( checkPointMsg,
111			"returned from the force calculation.\n" );
112			MPIcheckPoint();
113			#endif // is_mpi
114	mmeineke	597
115	mmeineke	377
116			}
117
118
119			void ForceFields::initFortran(int ljMixPolicy, int useReactionField ){
120
121			int isError;
122
123			isError = 0;
124			initFortranFF( &ljMixPolicy, &useReactionField, &isError );
125
126			if(isError){
127			sprintf( painCave.errMsg,
128			"ForceField error: There was an error initializing the forceField in fortran.\n" );
129			painCave.isFatal = 1;
130			simError();
131			}
132
133
134			#ifdef IS_MPI
135			sprintf( checkPointMsg, "ForceField successfully initialized the fortran component list.\n" );
136			MPIcheckPoint();
137			#endif // is_mpi
138
139			}