ParticlePairs2PotForceAdapter.h

 
#ifndef PARTICLEPAIRS2POTFORCEADAPTER_H_
#define PARTICLEPAIRS2POTFORCEADAPTER_H_
 
#include "molecules/potforce.h"
#include "particleContainer/handlerInterfaces/ParticlePairsHandler.h"
#include "io/RDF.h"
#include "Domain.h"
#include "WrapOpenMP.h"
 
class ParticlePairs2PotForceAdapter : public ParticlePairsHandler {
public:
    ParticlePairs2PotForceAdapter(Domain& domain) :
        _domain(domain), _virial(0.0), _upot6LJ(0.0), _upotXpoles(0.0), _myRF(0.0) {
        //this->_doRecordRDF = false;
 
        const int numThreads = mardyn_get_max_threads();
        Log::global_log->info() << "ParticlePairs2PotForceAdapter: allocate data for " << numThreads << " threads." << std::endl;
        _threadData.resize(numThreads);
        #if defined(_OPENMP)
        #pragma omp parallel
        #endif
        {
            PP2PFAThreadData * myown = new PP2PFAThreadData();
            const int myid = mardyn_get_thread_num();
            _threadData[myid] = myown;
        } // end pragma omp parallel
    }
 
    ~ParticlePairs2PotForceAdapter() {
        #if defined(_OPENMP)
        #pragma omp parallel
        #endif
        {
            const int myid = mardyn_get_thread_num();
            delete _threadData[myid];
        } // end pragma omp parallel
    }
 
    void init() {
        _virial = 0;
        _upot6LJ = 0;
        _upotXpoles = 0;
        _myRF = 0;
 
        #if defined(_OPENMP)
        #pragma omp parallel
        #endif
        {
            const int myid = mardyn_get_thread_num();
            _threadData[myid]->initComp2Param(_domain.getComp2Params());
            _threadData[myid]->clear();
        }
    }
 
    void finish() {
        double glob_virial = 0.0;
        double glob_upot6LJ = 0.0;
        double glob_upotXpoles = 0.0;
        double glob_myRF = 0.0;
 
        #if defined(_OPENMP)
        #pragma omp parallel reduction(+:glob_virial, glob_upot6LJ, glob_upotXpoles, glob_myRF)
        #endif
        {
            const int myid = mardyn_get_thread_num();
            glob_virial += _threadData[myid]->_virial;
            glob_upot6LJ += _threadData[myid]->_upot6LJ;
            glob_upotXpoles += _threadData[myid]->_upotXpoles;
            glob_myRF += _threadData[myid]->_myRF;
        } // end pragma omp parallel reduction
 
        _virial = glob_virial;
        _upot6LJ = glob_upot6LJ;
        _upotXpoles = glob_upotXpoles;
        _myRF = glob_myRF;
 
        _domain.setLocalUpot(_upot6LJ / 6. + _upotXpoles + _myRF);
        _domain.setLocalVirial(_virial + 3.0 * _myRF);
    }
 
    struct PP2PFAThreadData {
        PP2PFAThreadData() : _virial(0.0), _upot6LJ(0.0), _upotXpoles(0.0), _myRF(0.0), _comp2Param(0) {}
 
        ~~PP2PFAThreadData() {
            if(_comp2Param != 0) {
                delete _comp2Param;
                _comp2Param = 0;
            }
        }
 
        void initComp2Param(Comp2Param& c2p) {
            if (_comp2Param != 0) {
                delete _comp2Param;
            }
            _comp2Param = new Comp2Param(c2p);
        }
 
        void clear() {
            _virial = 0.0;
            _upot6LJ = 0.0;
            _upotXpoles = 0.0;
            _myRF = 0.0;
        }
 
        double _virial;
        double _upot6LJ;
        double _upotXpoles;
        double _myRF;
        Comp2Param * _comp2Param;
    };
 
    std::vector<PP2PFAThreadData *> _threadData;
 
    double processPair(Molecule& molecule1, Molecule& molecule2, double distanceVector[3], PairType pairType, double dd, bool calculateLJ = true) {
        const int tid = mardyn_get_thread_num();
        PP2PFAThreadData &my_threadData = *_threadData[tid];
 
        ParaStrm& params = (* my_threadData._comp2Param)(molecule1.componentid(), molecule2.componentid());
        params.reset_read();
 
        switch (pairType) {
 
            double dummy1, dummy2, dummy3, dummy4[3], Virial3[3];
            
            case MOLECULE_MOLECULE : 
//                if ( _rdf != NULL ) _rdf->observeRDF(molecule1, molecule2, dd); // moved to RDFCellProcessor
                PotForce(molecule1, molecule2, params, distanceVector, my_threadData._upot6LJ, my_threadData._upotXpoles, my_threadData._myRF, Virial3, calculateLJ );
                my_threadData._virial += 2*(Virial3[0]+Virial3[1]+Virial3[2]);
                return my_threadData._upot6LJ + my_threadData._upotXpoles;
            case MOLECULE_HALOMOLECULE : 
 
                PotForce(molecule1, molecule2, params, distanceVector, dummy1, dummy2, dummy3, dummy4, calculateLJ);
                return 0.0;
            case MOLECULE_MOLECULE_FLUID : 
                dummy1 = 0.0; // 6*U_LJ
                dummy2 = 0.0; // U_polarity
                dummy3 = 0.0; // U_dipole_reaction_field
 
                FluidPot(molecule1, molecule2, params, distanceVector, dummy1, dummy2, dummy3, calculateLJ);
                return dummy1 / 6.0 + dummy2 + dummy3;
            default:
                Simulation::exit(666);
        }
        return 0.0;
    }
 
 
//  void recordRDF() {
//      this->_doRecordRDF = true;
//  }
 
private:
    Domain& _domain;
 
    double _virial;
    double _upot6LJ;
    double _upotXpoles;
    double _myRF;
 
//  bool _doRecordRDF;
};
 
#endif /*PARTICLEPAIRS2POTFORCEADAPTER_H_*/