GBVIParameters.cpp

/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
 * Contributors: Pande Group
 *
 * 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 <math.h>
#include <sstream>
#include <string.h>

#include "openmm/OpenMMException.h"
#include "GBVIParameters.h"

using std::vector;
using namespace OpenMM;

/**---------------------------------------------------------------------------------------

    GBVIParameters constructor 

    @param numberOfAtoms       number of atoms

    --------------------------------------------------------------------------------------- */

GBVIParameters::GBVIParameters(int numberOfAtoms) : _numberOfAtoms(numberOfAtoms), 
                                                    _soluteDielectric(1.0),
                                                    _solventDielectric(78.3),
                                                    _electricConstant(-0.5*ONE_4PI_EPS0),
                                                    _cutoff(false),
                                                    _periodic(false),
                                                    _bornRadiusScalingMethod(0),
                                                    _quinticLowerLimitFactor(0.8),
                                                    _quinticUpperBornRadiusLimit(5.0) {

    _atomicRadii.resize(numberOfAtoms);
    _scaledRadii.resize(numberOfAtoms);
    _gammaParameters.resize(numberOfAtoms); 

}

/**---------------------------------------------------------------------------------------

    GBVIParameters destructor 

    --------------------------------------------------------------------------------------- */

GBVIParameters::~GBVIParameters() {
}

/**---------------------------------------------------------------------------------------

   Get number of atoms

   @return number of atoms

   --------------------------------------------------------------------------------------- */

int GBVIParameters::getNumberOfAtoms() const {
    return _numberOfAtoms;
}

/**---------------------------------------------------------------------------------------

   Get electric constant

   @return electric constant

   --------------------------------------------------------------------------------------- */

RealOpenMM GBVIParameters::getElectricConstant() const {
    return _electricConstant;
}

/**---------------------------------------------------------------------------------------

   Get solvent dielectric

   @return solvent dielectric

   --------------------------------------------------------------------------------------- */

RealOpenMM GBVIParameters::getSolventDielectric() const {
    return _solventDielectric;
}

/**---------------------------------------------------------------------------------------

   Set solvent dielectric

   @param solventDielectric solvent dielectric

   --------------------------------------------------------------------------------------- */

void GBVIParameters::setSolventDielectric(RealOpenMM solventDielectric) {
    _solventDielectric = solventDielectric;
}

/**---------------------------------------------------------------------------------------

   Get solute dielectric

   @return soluteDielectric

   --------------------------------------------------------------------------------------- */

RealOpenMM GBVIParameters::getSoluteDielectric() const {
    return _soluteDielectric;
}

/**---------------------------------------------------------------------------------------

   Set solute dielectric

   @param soluteDielectric solute dielectric

   --------------------------------------------------------------------------------------- */

void GBVIParameters::setSoluteDielectric(RealOpenMM soluteDielectric) {
    _soluteDielectric = soluteDielectric;
}

/**---------------------------------------------------------------------------------------

    Get AtomicRadii array

    @return array of atomic radii

    --------------------------------------------------------------------------------------- */

const vector<RealOpenMM>& GBVIParameters::getAtomicRadii() const {
    return _atomicRadii;
}

/**---------------------------------------------------------------------------------------

    Set AtomicRadii array

    @param atomicRadii vector of atomic radii

    --------------------------------------------------------------------------------------- */

void GBVIParameters::setAtomicRadii(const vector<RealOpenMM>& atomicRadii) {

    if (atomicRadii.size() == _atomicRadii.size()) {
        for (unsigned int ii = 0; ii < atomicRadii.size(); ii++) {
            _atomicRadii[ii] = atomicRadii[ii];
        }   
    } else {
        std::stringstream msg;
        msg << "GBVIParameters: input size for atomic radii does not agree w/ current size: input=";
        msg << atomicRadii.size();
        msg << " current size=" << _atomicRadii.size();
        throw OpenMM::OpenMMException(msg.str());
    }

}

/**---------------------------------------------------------------------------------------

    Return scaled radii
    @return array 

    --------------------------------------------------------------------------------------- */

const vector<RealOpenMM>& GBVIParameters::getScaledRadii() const {
    return _scaledRadii;
}

/**---------------------------------------------------------------------------------------

    Set scaled radii

    @param scaledRadii  scaledRadii

    --------------------------------------------------------------------------------------- */

void GBVIParameters::setScaledRadii(const vector<RealOpenMM>& scaledRadii) {

    if (scaledRadii.size() == _scaledRadii.size()) {
        for (unsigned int ii = 0; ii < scaledRadii.size(); ii++) {
            _scaledRadii[ii] = scaledRadii[ii];
        }
    } else {
        std::stringstream msg;
        msg << "GBVIParameters: input size for scaled radii does not agree w/ current size: input=";
        msg << scaledRadii.size();
        msg << " current size=" << _scaledRadii.size();
        throw OpenMM::OpenMMException(msg.str());
    }

}

/**---------------------------------------------------------------------------------------

    Return gamma parameters
    If not previously set, allocate space

    @return array 

    --------------------------------------------------------------------------------------- */

const vector<RealOpenMM>& GBVIParameters::getGammaParameters() const {
    return _gammaParameters;
}

/**---------------------------------------------------------------------------------------

    Set gamma parameters

    @param gammas  gammas

    --------------------------------------------------------------------------------------- */

void GBVIParameters::setGammaParameters(const vector<RealOpenMM>& gammas) {

    if (gammas.size() == _gammaParameters.size()) {
        for (unsigned int ii = 0; ii < gammas.size(); ii++) {
            _gammaParameters[ii] = gammas[ii];
        }
    } else {
        std::stringstream msg;
        msg << "GBVIParameters: input size for gammas does not agree w/ current size: input=";
        msg << gammas.size();
        msg << " current size=" << _gammaParameters.size();
        throw OpenMM::OpenMMException(msg.str());
    }

}

/**---------------------------------------------------------------------------------------

      Set the force to use a cutoff.

      @param distance            the cutoff distance

      --------------------------------------------------------------------------------------- */

void GBVIParameters::setUseCutoff(RealOpenMM distance) {

     _cutoff          = true;
     _cutoffDistance = distance;
}

/**---------------------------------------------------------------------------------------

      Get whether to use a cutoff.

      --------------------------------------------------------------------------------------- */

bool GBVIParameters::getUseCutoff() {
     return _cutoff;
}

/**---------------------------------------------------------------------------------------

      Get the cutoff distance.

      --------------------------------------------------------------------------------------- */

RealOpenMM GBVIParameters::getCutoffDistance() {
     return _cutoffDistance;
}

/**---------------------------------------------------------------------------------------

      Set the force to use periodic boundary conditions.  This requires that a cutoff has
      also been set, and the smallest side of the periodic box is at least twice the cutoff
      distance.

      @param vectors    the vectors defining the periodic box

      --------------------------------------------------------------------------------------- */

void GBVIParameters::setPeriodic(RealVec* vectors) {

    assert(_cutoff);
    assert(vectors[0][0] >= 2.0*_cutoffDistance);
    assert(vectors[1][1] >= 2.0*_cutoffDistance);
    assert(vectors[2][2] >= 2.0*_cutoffDistance);
    _periodic = true;
    _periodicBoxVectors[0] = vectors[0];
    _periodicBoxVectors[1] = vectors[1];
    _periodicBoxVectors[2] = vectors[2];
}

/**---------------------------------------------------------------------------------------

      Get whether to use periodic boundary conditions.

      --------------------------------------------------------------------------------------- */

bool GBVIParameters::getPeriodic() {
     return _periodic;
}

/**---------------------------------------------------------------------------------------

      Get the periodic box dimension

      --------------------------------------------------------------------------------------- */

const OpenMM::RealVec* GBVIParameters::getPeriodicBox() {
     return _periodicBoxVectors;
}

/**---------------------------------------------------------------------------------------

    Get tau prefactor

    @return (1/e1 - 1/e0), where e1 = solute dielectric, e0 = solvent dielectric

    --------------------------------------------------------------------------------------- */

RealOpenMM GBVIParameters::getTau() const {

    // ---------------------------------------------------------------------------------------

    static const RealOpenMM zero = 0.0;
    static const RealOpenMM one  = 1.0;

    // ---------------------------------------------------------------------------------------

    RealOpenMM tau;
    if (getSoluteDielectric() != zero && getSolventDielectric() != zero) {
        tau = (one/getSoluteDielectric()) - (one/getSolventDielectric());
    } else {
        tau = zero;
    }   

    return tau;
}

/**---------------------------------------------------------------------------------------

    Get bornRadiusScalingMethod

    @return bornRadiusScalingMethod

    --------------------------------------------------------------------------------------- */

int GBVIParameters::getBornRadiusScalingMethod() const {
    return _bornRadiusScalingMethod;
}

/**---------------------------------------------------------------------------------------

    Set bornRadiusScalingMethod 

    @param bornRadiusScalingMethod bornRadiusScalingMethod

    --------------------------------------------------------------------------------------- */

void GBVIParameters::setBornRadiusScalingMethod(int bornRadiusScalingMethod) {
    _bornRadiusScalingMethod    = bornRadiusScalingMethod;
}

/**---------------------------------------------------------------------------------------

    Get quinticLowerLimitFactor

    @return quinticLowerLimitFactor

    --------------------------------------------------------------------------------------- */

RealOpenMM GBVIParameters::getQuinticLowerLimitFactor() const {
    return _quinticLowerLimitFactor;
}

/**---------------------------------------------------------------------------------------

    Set quinticLowerLimitFactor 

    @param quinticLowerLimitFactor quinticLowerLimitFactor

    --------------------------------------------------------------------------------------- */

void GBVIParameters::setQuinticLowerLimitFactor(RealOpenMM quinticLowerLimitFactor) {
    _quinticLowerLimitFactor    = quinticLowerLimitFactor;
}

/**---------------------------------------------------------------------------------------

    Get quinticUpperBornRadiusLimit

    @return quinticUpperBornRadiusLimit

    --------------------------------------------------------------------------------------- */

RealOpenMM GBVIParameters::getQuinticUpperBornRadiusLimit() const {
    return _quinticUpperBornRadiusLimit;
}

/**---------------------------------------------------------------------------------------

    Set quinticUpperBornRadiusLimit 

    @param quinticUpperBornRadiusLimit quinticUpperBornRadiusLimit

    --------------------------------------------------------------------------------------- */

void GBVIParameters::setQuinticUpperBornRadiusLimit(RealOpenMM quinticUpperBornRadiusLimit) {
    _quinticUpperBornRadiusLimit    = quinticUpperBornRadiusLimit;
}