tfPotential.h

Go to the documentation of this file.

/*******************************************************************************
 * This file is part of mdcore.
 * Coypright (c) 2010 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
 * Coypright (c) 2017 Andy Somogyi (somogyie at indiana dot edu)
 * Copyright (c) 2022-2024 T.J. Sego
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

 
#ifndef _MDCORE_INCLUDE_TFPOTENTIAL_H_
#define _MDCORE_INCLUDE_TFPOTENTIAL_H_
 
#include "tf_platform.h"
#include "tf_fptype.h"
#include <io/tf_io.h>
 
#include <limits>
#include <utility>
#include <vector>
 
/* some constants */
#define potential_degree                    5
#define potential_chunk                     (potential_degree+3)
#define potential_ivalsa                    1
#define potential_ivalsb                    10
#define potential_N                         100
#define potential_align                     64
#define potential_ivalsmax                  3048
 
#define potential_escale                    (0.079577471545947667882)
// #define potential_escale                    1.0
 
 
namespace TissueForge {
 
 
    /* potential flags */
 
    enum PotentialFlags {
        POTENTIAL_NONE            = 0,
        POTENTIAL_LJ126           = 1 << 0,
        POTENTIAL_EWALD           = 1 << 1,
        POTENTIAL_COULOMB         = 1 << 2,
        POTENTIAL_SINGLE          = 1 << 3,

        POTENTIAL_R2              = 1 << 4,

        POTENTIAL_R               = 1 << 5,

        POTENTIAL_ANGLE           = 1 << 6,

        POTENTIAL_HARMONIC        = 1 << 7,
 
        POTENTIAL_DIHEDRAL        = 1 << 8,

        POTENTIAL_SWITCH          = 1 << 9,
 
        POTENTIAL_REACTIVE        = 1 << 10,

        POTENTIAL_SCALED          = 1 << 11,

        POTENTIAL_SHIFTED         = 1 << 12,

        POTENTIAL_BOUND           = 1 << 13,
 
        // sum of constituent potentials
        POTENTIAL_SUM             = 1 << 14, 

        POTENTIAL_PERIODIC        = 1 << 15, 
 
        /* Coulomb reciprocated */
        POTENTIAL_COULOMBR        = 1 << 16
    };
 
    enum PotentialKind {
        // standard interpolated potential kind
        POTENTIAL_KIND_POTENTIAL,
        
        // dissipative particle dynamics kind
        POTENTIAL_KIND_DPD, 
 
        // explicit potential by particles
        POTENTIAL_KIND_BYPARTICLES,
 
        // combination of two constituent potentials
        POTENTIAL_KIND_COMBINATION
    };
 

    typedef void (*PotentialEval_ByParticle) (
        struct Potential *p, 
        struct Particle *part_i, 
        FPTYPE *dx, 
        FPTYPE r2, 
        FPTYPE *e, 
        FPTYPE *f
    );

    typedef void (*PotentialEval_ByParticles) (
        struct Potential *p, 
        struct Particle *part_i, 
        struct Particle *part_j, 
        FPTYPE *dx, 
        FPTYPE r2, 
        FPTYPE *e, 
        FPTYPE *f
    );

    typedef void (*PotentialEval_ByParticles3)(
        struct Potential *p, 
        struct Particle *part_i, 
        struct Particle *part_j, 
        struct Particle *part_k, 
        FPTYPE ctheta, 
        FPTYPE *e, 
        FPTYPE *fi, 
        FPTYPE *fk
    );

    typedef void (*PotentialEval_ByParticles4)(
        struct Potential *p, 
        struct Particle *part_i, 
        struct Particle *part_j, 
        struct Particle *part_k, 
        struct Particle *part_l, 
        FPTYPE cphi, 
        FPTYPE *e, 
        FPTYPE *fi, 
        FPTYPE *fl
    );
 
    typedef struct Potential* (*PotentialCreate) (
        struct Potential *partial_potential,
        struct ParticleType *a, 
        struct ParticleType *b
    );

    typedef void (*PotentialClear) (struct Potential* p);
 

    struct CAPI_EXPORT Potential {
        uint32_t kind;

        uint32_t flags;
 

        FPTYPE alpha[4];

        FPTYPE *c;
 
        FPTYPE r0_plusone;

        FPTYPE a, b;

        FPTYPE mu;

        FPTYPE offset[3];

        int n;
 
        PotentialCreate create_func;
        PotentialClear clear_func;
 
        PotentialEval_ByParticle eval_bypart;
        PotentialEval_ByParticles eval_byparts;
        PotentialEval_ByParticles3 eval_byparts3;
        PotentialEval_ByParticles4 eval_byparts4;
 
        Potential *pca, *pcb;

        const char* name;
 
        Potential();
 
        FPTYPE operator()(const FPTYPE &r, const FPTYPE &r0=-1.0);
        FPTYPE operator()(const std::vector<FPTYPE>& r);
        FPTYPE operator()(struct ParticleHandle* pi, const FVector3 &pt);
        FPTYPE operator()(struct ParticleHandle* pi, struct ParticleHandle* pj);
        FPTYPE operator()(struct ParticleHandle* pi, struct ParticleHandle* pj, struct ParticleHandle* pk);
        FPTYPE operator()(struct ParticleHandle* pi, struct ParticleHandle* pj, struct ParticleHandle* pk, struct ParticleHandle* pl);
        FPTYPE force(FPTYPE r, FPTYPE ri=-1.0, FPTYPE rj=-1.0);
        std::vector<FPTYPE> force(const std::vector<FPTYPE>& r);
        std::vector<FPTYPE> force(struct ParticleHandle* pi, const FVector3 &pt);
        std::vector<FPTYPE> force(struct ParticleHandle* pi, struct ParticleHandle* pj);
        std::pair<std::vector<FPTYPE>, std::vector<FPTYPE> > force(struct ParticleHandle* pi, struct ParticleHandle* pj, struct ParticleHandle* pk);
        std::pair<std::vector<FPTYPE>, std::vector<FPTYPE> > force(struct ParticleHandle* pi, struct ParticleHandle* pj, struct ParticleHandle* pk, struct ParticleHandle* pl);
 
        std::vector<Potential*> constituents();
 
        Potential& operator+(const Potential& rhs);

        virtual std::string toString();

        static Potential *fromString(const std::string &str);

        static Potential *lennard_jones_12_6(FPTYPE min, FPTYPE max, FPTYPE A, FPTYPE B, FPTYPE *tol=NULL);

        static Potential *lennard_jones_12_6_coulomb(FPTYPE min, FPTYPE max, FPTYPE A, FPTYPE B, FPTYPE q, FPTYPE *tol=NULL);

        static Potential *ewald(FPTYPE min, FPTYPE max, FPTYPE q, FPTYPE kappa, FPTYPE *tol=NULL, unsigned int *periodicOrder=NULL);

        static Potential *coulomb(FPTYPE q, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL, unsigned int *periodicOrder=NULL);

        static Potential* coulombR(FPTYPE q, FPTYPE kappa, FPTYPE min, FPTYPE max, unsigned int* modes=NULL);

        static Potential *harmonic(FPTYPE k, FPTYPE r0, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL);

        static Potential *linear(FPTYPE k, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL);

        static Potential *harmonic_angle(FPTYPE k, FPTYPE theta0, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL);

        static Potential *harmonic_dihedral(FPTYPE k, FPTYPE delta, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL);

        static Potential *cosine_dihedral(FPTYPE k, int n, FPTYPE delta, FPTYPE *tol=NULL);

        static Potential *well(FPTYPE k, FPTYPE n, FPTYPE r0, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL);

        static Potential *glj(FPTYPE e, FPTYPE *m=NULL, FPTYPE *n=NULL, FPTYPE *k=NULL, FPTYPE *r0=NULL, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL, bool *shifted=NULL);

        static Potential *morse(FPTYPE *d=NULL, FPTYPE *a=NULL, FPTYPE *r0=NULL, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL, bool *shifted=NULL);

        static Potential *overlapping_sphere(FPTYPE *mu=NULL, FPTYPE *kc=NULL, FPTYPE *kh=NULL, FPTYPE *r0=NULL, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL);
        
        static Potential *power(FPTYPE *k=NULL, FPTYPE *r0=NULL, FPTYPE *alpha=NULL, FPTYPE *min=NULL, FPTYPE *max=NULL, FPTYPE *tol=NULL);

        static Potential *dpd(FPTYPE *alpha=NULL, FPTYPE *gamma=NULL, FPTYPE *sigma=NULL, FPTYPE *cutoff=NULL, bool *shifted=NULL);

        static Potential *custom(
            FPTYPE min, 
            FPTYPE max, 
            FPTYPE (*f)(FPTYPE), 
            FPTYPE (*fp)(FPTYPE), 
            FPTYPE (*f6p)(FPTYPE), 
            FPTYPE *tol=NULL, 
            uint32_t *flags=NULL
        );
 
        FPTYPE getMin();
        FPTYPE getMax();
        FPTYPE getCutoff();
        std::pair<FPTYPE, FPTYPE> getDomain();
        int getIntervals();
        bool getBound();
        void setBound(const bool &_bound);
        FPTYPE getR0();
        void setR0(const FPTYPE &_r0);
        bool getShifted();
        bool getPeriodic();
        bool getRSquare();
 
    };
 

    CAPI_DATA(struct Potential) potential_null;
 
 
    /* associated functions */

    CAPI_FUNC(void) potential_clear(struct Potential *p);

    CAPI_FUNC(HRESULT) potential_init(
        struct Potential *p, 
        FPTYPE (*f)(FPTYPE),
        FPTYPE (*fp)(FPTYPE), 
        FPTYPE (*f6p)(FPTYPE),
        FPTYPE a, 
        FPTYPE b, 
        FPTYPE tol
    );
 
    CAPI_FUNC(HRESULT) potential_getcoeffs(
        FPTYPE (*f)(FPTYPE), 
        FPTYPE (*fp)(FPTYPE),
        FPTYPE *xi, 
        int n, 
        FPTYPE *c, 
        FPTYPE *err
    );
 
    CAPI_FUNC(FPTYPE) potential_getalpha(FPTYPE (*f6p)(FPTYPE), FPTYPE a, FPTYPE b);
 
};
 
#endif // _MDCORE_INCLUDE_TFPOTENTIAL_H_