ForceFields

ForceField.hpp

#ifndef FORCEFIELD_HPP
#define FORCEFIELD_HPP
#include "vector.hpp"
#include "atom.hpp"
#include "system.hpp"

// I moved all code to one file for gist easier sharing ;)

// change to potential energy??
class ForceField {
    
    public:
        // force after "interaction"
        virtual Vector3d calc_force(Atom* atom) = 0;
        // potential energy
        virtual double calc_u(Atom* atom) = 0;
        // representation of force_field if any
        virtual string get_representation();
        string const name;
        
        ForceField(string name):
            name(name) {}
};

string ForceField::get_representation() {
    return "{}";
}

class SoftWall : public ForceField
{
    private:
        Hyperplane<double, 3> *plane;
        
        double width;
        double softness; // hardness essentially(!)
        
    public:
        SoftWall(Vector3d normal, Vector3d position):
            SoftWall(normal, position, 1, 1) {}
        
        SoftWall(Vector3d normal, double distance):
            SoftWall(normal, distance, 1, 1) {}
            
        // position: a position of a point which should belong to the wall
        SoftWall(Vector3d normal, Vector3d position, double width, double softness):
            ForceField("soft_wall"), width(width), softness(softness)
            {
                // "Construct a plane from its normal n and a point e onto the plane."
                plane = new Hyperplane<double, 3>(normal, position);
            }
            
        // distance: distance from origin centre
        SoftWall(Vector3d normal, double distance, double width, double softness):
            ForceField("soft_wall"), width(width), softness(softness)
            {
                // "Constructs a plane from its normal n and distance to the origin d
                // such that the algebraic equation of the plane is $n \cdot x + d = 0$."
                plane = new Hyperplane<double, 3>(normal, distance);
            }
            
        double calc_u(Atom* atom) override;
        Vector3d calc_force(Atom* atom) override;
        string get_representation() override;
};


double SoftWall::calc_u(Atom* atom)
{
    double distance = plane->absDistance(atom->position);
    
    if (distance <= width)
    {
        // 1/2 * (w-sqrt(x^2+y^2+z^2))^2*s
        return 0.5 * pow(width - distance, 2) * softness;
    }
    return 0;
};

#include "format.hpp"
string SoftWall::get_representation()
{
    using namespace fmt;

    auto normal = this->plane->normal();

    return format(
        "{{"
            "\"type\": \"wall\","
            " \"width\": {width},"
            " \"x\": {x},"
            " \"y\": {y},"
            " \"z\": {z},"
            " \"distance\": {distance}"
        "}}",
        arg("width", this->width),
        arg("x", normal.x()),
        arg("y", normal.y()),
        arg("z", normal.z()),
        arg("distance", this->plane->offset())
    );
}

Vector3d SoftWall::calc_force(Atom* atom)
{
    Vector3d closest_on_plane = plane->projection(atom->position);

    Vector3d projection = closest_on_plane - atom->position;
    
    double distance = projection.norm();
    
    // assert distance == plane->absDistance(atom.position);
    
    if (distance <= width)
    {
        //  return - (2 * projection * softness);
        // derivative: $\vec{Force}(\vec{r}) = - \nabla \text{Potential energy}_{\vec{r}} $
        // derivative of: 1/2 * (w-sqrt(x^2+y^2+z^2))^2*s
        // sx(√x2+z2+y2−w)/√x2+z2+y2
        return softness * projection * (distance - width) / distance;
    }
    return Zero3d;
};

/*
 * Simplified Lennard-Jones potential,
 * @see AB form on https://en.wikipedia.org/wiki/Lennard-Jones_potential#AB_form
 * @see https://chem.libretexts.org/Core/Physical_and_Theoretical_Chemistry/Physical_Properties_of_Matter/Atomic_and_Molecular_Properties/Intermolecular_Forces/Specific_Interactions/Lennard-Jones_Potential
 * DOES NOT PLAY WELL
 */
class VanDerWaals: public ForceField
{
    private:
       System* system;
    public:
        VanDerWaals(System* system):
            ForceField("van_der_waals"), system(system) {}
        
        double calc_u(Atom* atom) override;
        Vector3d calc_force(Atom* atom) override;
};

double VanDerWaals::calc_u(Atom* atom)
{
    double energy = 0;
    
    for(Atom* other_atom : this->system->atoms)
    {
        if(other_atom == atom)
            continue;
        
        // how close two atom can get to each other (A)
        double sigma = atom->radius() + other_atom->radius();

        Vector3d difference = atom->position - other_atom->position;
        double distance = difference.norm();

        // cut-off; sometimes 12 * sigma. Should be somehow smoothed.
        if(distance > 6 * sigma)
            continue;

        // how strongly two atoms attract each other (kj/mol)
        double epsilon = 1;

        // (Pauli) repulsion
        double A = pow(sigma / distance, 12);

        // (London dispersion) attraction
        double B = pow(sigma / distance, 6);
        
        energy += 4 * epsilon * (A - B);
    }
    return energy;
    
}

Vector3d VanDerWaals::calc_force(Atom* atom)
{
    
    /* Force equals to $F = - \nabla U(r)$
     * Formula for U is given and explained in calc_u. Partial derivative with respect to x equals:
     * $\frac{\partial \:}{\partial \:x}(\frac{A}{(\sqrt{x^2+y^2+z^2})^{12}}) = \frac{12Ax}{(x^2+y^2+z^2)^7}$
     * for the attractive component and
     * $\frac{\partial \:}{\partial \:x}(\frac{B}{(\sqrt{x^2+y^2+z^2})^{7}}) = \frac{7Bx}{\sqrt{(x^2+y^2+z^2)}^{9}$
     */
    Vector3d force = Zero3d;
    
    for(Atom* other_atom : this->system->atoms)
    {
        if(other_atom == atom)
            continue;
        
        // how close two atom can get to each other (A)
        double sigma = atom->radius() + other_atom->radius();
        
        Vector3d difference = atom->position - other_atom->position;
        double distance = difference.norm();
        
        // cut-off; sometimes 12 * sigma. Should be somehow smoothed.
        if(distance > 6 * sigma)
            continue;
        
        // how strongly two atoms attract each other (kj/mol)
        double epsilon = 1;
        
        force += difference * 48 * epsilon / pow(sigma, 2) * (pow(sigma / distance, 14) - 0.5 * pow(sigma / distance, 8));
    }
    return force;
}



class Coulomb: public ForceField
{
    private:
       System* system;
    public:
        Coulomb(System* system):
            ForceField("coulomb"), system(system) {}
        
        double calc_u(Atom* atom) override;
        Vector3d calc_force(Atom* atom) override;
};

double Coulomb::calc_u(Atom* atom)
{
    double energy = 0;
    
    for(Atom* other_atom : this->system->atoms)
    {
        if(other_atom == atom)
            continue;
        
        Vector3d difference = atom->position - other_atom->position;
        double distance = difference.norm();
        energy += atom->charge() * other_atom->charge() / distance;
    }
    return energy;    
}

Vector3d Coulomb::calc_force(Atom *atom)
{
    /*
     * To check that sum of partial derivatives:
     *      $\frac{q_1 q_2 x}{\sqrt{x^2 + y^2 +z^2}^3} + \partial_y + partial_z$
     * equals force from Coulumb Law:
     *   raise each fraction (each partial) to power of two:
     *      $\frac{q_1^2 q_2^2 x^2}{\sqrt{x^2 + y^2 +z^2}^6} + \partial_y^2 + partial_z^2$
     *   add the fractions, move $q_1 q_2$ ahead of parens:
     *      $\frac{(q_1^2 q_2^2)(x^2 + y^2 + z^2)}{\sqrt{x^2 + y^2 +z^2}^6}$
     *   change sum of squares to square root:
     *      $\frac{(q_1 q_2)^2\sqrt{x^2 + y^2 + z^2}^2}{\sqrt{x^2 + y^2 +z^2}^6}$
     *   reduce:
     *      $\frac{(q_1 q_2)^2}{\sqrt{x^2 + y^2 +z^2}^4}$
     *   and again simplify too:
     *      $\frac{q_1 q_2}{\sqrt{x^2 + y^2 +z^2}^2}$
     * which is exactly formulation of Coulumb's Law
     */
    
    Vector3d force = Zero3d;
    
    for(Atom* other_atom : this->system->atoms)
    {
        if(other_atom == atom)
            continue;
        
        Vector3d difference = atom->position - other_atom->position;
        double distance = difference.norm();
        double qq = atom->charge() * other_atom->charge();
        force += atom->position * qq / pow(distance, 3);
    }
    
    return force;
}
#endif