AmoebaMultipoleForce.h

#ifndef OPENMM_AMOEBA_MULTIPOLE_FORCE_H_
#define OPENMM_AMOEBA_MULTIPOLE_FORCE_H_

/* -------------------------------------------------------------------------- *
 *                              OpenMMAmoeba                                  *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
 * Authors: Mark Friedrichs, Peter Eastman                                    *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include "openmm/Force.h"
#include "openmm/OpenMMException.h"
#include "internal/windowsExportAmoeba.h"
#include "openmm/Vec3.h"

#include <sstream>
#include <vector>

namespace OpenMM {

class OPENMM_EXPORT_AMOEBA AmoebaMultipoleForce : public Force {

public:

    enum NonbondedMethod {

        NoCutoff = 0,

        PME = 1
    };

    enum PolarizationType {

        Mutual = 0,

        Direct = 1
    };

    enum MultipoleAxisTypes { ZThenX = 0, Bisector = 1, ZBisect = 2, ThreeFold = 3, ZOnly = 4, NoAxisType = 5, LastAxisTypeIndex = 6 };

    enum CovalentType {
                          Covalent12 = 0, Covalent13 = 1, Covalent14 = 2, Covalent15 = 3,
                          PolarizationCovalent11 = 4, PolarizationCovalent12 = 5, PolarizationCovalent13 = 6, PolarizationCovalent14 = 7, CovalentEnd = 8 };

    AmoebaMultipoleForce();

    int getNumMultipoles() const {
        return multipoles.size();
    }

    NonbondedMethod getNonbondedMethod() const;

    void setNonbondedMethod(NonbondedMethod method);

    PolarizationType getPolarizationType() const;

    void setPolarizationType(PolarizationType type);

    double getCutoffDistance() const;

    void setCutoffDistance(double distance);

    double getAEwald() const;

    void setAEwald(double aewald);

    int getPmeBSplineOrder() const;

    void getPmeGridDimensions(std::vector<int>& gridDimension) const;

    void setPmeGridDimensions(const std::vector<int>& gridDimension);

    int addMultipole(double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, int axisType,
                     int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity);

    void getMultipoleParameters(int index, double& charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole,
                                int& axisType, int& multipoleAtomZ, int& multipoleAtomX, int& multipoleAtomY, double& thole, double& dampingFactor, double& polarity) const;

    void setMultipoleParameters(int index, double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole,
                                int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity);

    void setCovalentMap(int index, CovalentType typeId, const std::vector<int>& covalentAtoms);

    void getCovalentMap(int index, CovalentType typeId, std::vector<int>& covalentAtoms) const;

    void getCovalentMaps(int index, std::vector < std::vector<int> >& covalentLists) const;

    int getMutualInducedMaxIterations(void) const;

    void setMutualInducedMaxIterations(int inputMutualInducedMaxIterations);

    double getMutualInducedTargetEpsilon(void) const;

    void setMutualInducedTargetEpsilon(double inputMutualInducedTargetEpsilon);

    double getEwaldErrorTolerance() const;
    void setEwaldErrorTolerance(double tol);

    void getInducedDipoles(Context& context, std::vector<Vec3>& dipoles);
    
    void getElectrostaticPotential(const std::vector< Vec3 >& inputGrid,
                                    Context& context, std::vector< double >& outputElectrostaticPotential);

    void getSystemMultipoleMoments(Context& context, std::vector< double >& outputMultipoleMoments);
    void updateParametersInContext(Context& context);
    bool usesPeriodicBoundaryConditions() const {
        return nonbondedMethod == AmoebaMultipoleForce::PME;
    }
protected:
    ForceImpl* createImpl() const;
private:
    NonbondedMethod nonbondedMethod;
    PolarizationType polarizationType;
    double cutoffDistance;
    double aewald;
    int pmeBSplineOrder;
    std::vector<int> pmeGridDimension;
    int mutualInducedMaxIterations;
    double mutualInducedTargetEpsilon;
    double scalingDistanceCutoff;
    double electricConstant;
    double ewaldErrorTol;
    class MultipoleInfo;
    std::vector<MultipoleInfo> multipoles;
};

class AmoebaMultipoleForce::MultipoleInfo {
public:

    int axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
    double charge, thole, dampingFactor, polarity;

    std::vector<double> molecularDipole;
    std::vector<double> molecularQuadrupole;
    std::vector< std::vector<int> > covalentInfo;

    MultipoleInfo() {
        axisType = multipoleAtomZ = multipoleAtomX = multipoleAtomY = -1;
        charge   = thole          = dampingFactor  = 0.0;

        molecularDipole.resize(3);
        molecularQuadrupole.resize(9);

    }

    MultipoleInfo(double charge, const std::vector<double>& inputMolecularDipole, const std::vector<double>& inputMolecularQuadrupole,
                   int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) :
        charge(charge), axisType(axisType), multipoleAtomZ(multipoleAtomZ), multipoleAtomX(multipoleAtomX), multipoleAtomY(multipoleAtomY),
        thole(thole), dampingFactor(dampingFactor), polarity(polarity) {

       covalentInfo.resize(CovalentEnd);

       molecularDipole.resize(3);
       molecularDipole[0]          = inputMolecularDipole[0];
       molecularDipole[1]          = inputMolecularDipole[1];
       molecularDipole[2]          = inputMolecularDipole[2];

       molecularQuadrupole.resize(9);
       molecularQuadrupole[0]      = inputMolecularQuadrupole[0];
       molecularQuadrupole[1]      = inputMolecularQuadrupole[1];
       molecularQuadrupole[2]      = inputMolecularQuadrupole[2];
       molecularQuadrupole[3]      = inputMolecularQuadrupole[3];
       molecularQuadrupole[4]      = inputMolecularQuadrupole[4];
       molecularQuadrupole[5]      = inputMolecularQuadrupole[5];
       molecularQuadrupole[6]      = inputMolecularQuadrupole[6];
       molecularQuadrupole[7]      = inputMolecularQuadrupole[7];
       molecularQuadrupole[8]      = inputMolecularQuadrupole[8];
    }
};

} // namespace OpenMM

#endif /*OPENMM_AMOEBA_MULTIPOLE_FORCE_H_*/