etomica.potential.P2HSDipole.java Source code

Java tutorial

Introduction

Here is the source code for etomica.potential.P2HSDipole.java

Source

/* This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

package etomica.potential;

import org.apache.commons.math3.special.Erf;

import etomica.api.IAtom;
import etomica.api.IAtomList;
import etomica.api.IBoundary;
import etomica.api.IBox;
import etomica.api.IMolecule;
import etomica.api.IMoleculeList;
import etomica.api.IVector;
import etomica.api.IVectorMutable;
import etomica.atom.IAtomOriented;
import etomica.exception.MethodNotImplementedException;
import etomica.space.ISpace;
import etomica.space.Tensor;

/**
 * Hard sphere molecule with a dipole sitting at the center.
 *
 * @author Weisong & Shu
 * date:May 2015
 * 
 */
public class P2HSDipole extends PotentialMolecular implements IPotentialMolecularSecondDerivative {
    //public class P2HSDipole extends PotentialMolecular implements IPotentialMolecularTorque  {
    public P2HSDipole(ISpace space) {
        this(space, 1, 1);
    }

    public P2HSDipole(ISpace space, double dipole) {
        this(space, 1, dipole);
    }

    public P2HSDipole(ISpace space, double sigma, double dipole) {
        this(space, 1, dipole, Double.POSITIVE_INFINITY);
    }

    public P2HSDipole(ISpace space, double sigma, double dipole, double rCut) {
        super(2, space);
        setSigma(sigma);
        gradient = new IVectorMutable[2];
        gradient[0] = space.makeVector();
        gradient[1] = space.makeVector();
        torque = new IVectorMutable[2];
        torque[0] = space.makeVector();
        torque[1] = space.makeVector();
        secondDerivative = new Tensor[3];
        this.secondDerivative[0] = space.makeTensor();
        this.secondDerivative[1] = space.makeTensor();
        this.secondDerivative[2] = space.makeTensor();
        a = new IVectorMutable[3];
        a[0] = space.makeVector();
        a[1] = space.makeVector();
        a[2] = space.makeVector();

        dr = space.makeVector();
        drunit = space.makeVector();
        runit = space.makeVector();
        work = space.makeVector();
        gradientAndTorque = new IVectorMutable[][] { gradient, torque };

        setDipole(dipole);
        this.rCut = rCut;
    }

    public void setBox(IBox box) {
        boundary = box.getBoundary();
    }

    public double getRange() {
        return Double.POSITIVE_INFINITY;
    }

    public double energy(IMoleculeList pair) {
        IMolecule molecule1 = pair.getMolecule(0);
        IMolecule molecule2 = pair.getMolecule(1);
        IAtomList atomList1 = molecule1.getChildList();
        IAtomList atomList2 = molecule2.getChildList();

        IAtomOriented atom1 = (IAtomOriented) atomList1.getAtom(0);
        IAtomOriented atom2 = (IAtomOriented) atomList2.getAtom(0);
        // dr is a r1-r2
        dr.Ev1Mv2(atom1.getPosition(), atom2.getPosition());
        boundary.nearestImage(dr);
        double r2 = dr.squared();

        if (r2 > rCut * rCut) {
            return 0;
        }

        if (r2 < sigma2) { // hard core
            return Double.POSITIVE_INFINITY;
        }
        // normalize dr, the vector between the molecules
        dr.normalize();

        // v1 (unit vector) is the orientation of molecule 1: dipole1 direction
        IVector v1 = atom1.getOrientation().getDirection();
        // v2 (unit vector) is the orientation of molecule 2: dipole1 direction
        IVector v2 = atom2.getOrientation().getDirection();

        // cos(dipole 1 and dipole 2)=cos(v1 and v2)
        double cos_D1_D2 = v1.dot(v2);
        //cos(dipole 1 and r12)
        double cos_D1_r = v1.dot(dr);
        //cos(r12 and dipole 2)
        double cos_r_D2 = dr.dot(v2);
        double r12Magnitude = Math.sqrt(r2);
        double ener = dipole * dipole * (cos_D1_D2 - 3.0 * cos_D1_r * cos_r_D2);
        ener = ener / r12Magnitude / r2;
        return ener;
    }

    public double getSigma() {
        return sigma;
    }

    public final void setSigma(double s) {
        sigma = s;
        sigma2 = s * s;
    }

    public void setDipole(double moment) {
        dipole = moment;
    }

    public IVector[] gradient(IMoleculeList pair, Tensor pressureTensor) {
        return gradient(pair);
    }

    public IVector[] gradient(IMoleculeList pair) {
        // do extra work to calculate torque
        gradientAndTorque(pair);
        return gradient;
    }

    public double virial(IMoleculeList atoms) {
        gradient(atoms);
        IMolecule molecule1 = atoms.getMolecule(0);
        IMolecule molecule2 = atoms.getMolecule(1);
        IAtomList atomList1 = molecule1.getChildList();
        IAtomList atomList2 = molecule2.getChildList();

        IAtomOriented atom1 = (IAtomOriented) atomList1.getAtom(0);
        IAtomOriented atom2 = (IAtomOriented) atomList2.getAtom(0);

        // LJ contributation

        dr.Ev1Mv2(atom2.getPosition(), atom1.getPosition());
        boundary.nearestImage(dr);
        double v = gradient[1].dot(dr);
        if (Double.isInfinite(v)) {
            throw new RuntimeException("oops " + v);
        }
        return gradient[1].dot(dr);
    }

    public IVector[][] gradientAndTorque(IMoleculeList atoms) {
        IMolecule molecule1 = atoms.getMolecule(0);
        IMolecule molecule2 = atoms.getMolecule(1);
        IAtomList atomList1 = molecule1.getChildList();
        IAtomList atomList2 = molecule2.getChildList();

        IAtomOriented atom1 = (IAtomOriented) atomList1.getAtom(0);
        IAtomOriented atom2 = (IAtomOriented) atomList2.getAtom(0);

        dr.Ev1Mv2(atom2.getPosition(), atom1.getPosition());
        boundary.nearestImage(dr);
        double r2 = dr.squared();
        if (r2 > rCut * rCut) {
            gradient[0].E(0);
            gradient[1].E(0);
            torque[0].E(0);
            torque[1].E(0);
            return gradientAndTorque;
        }

        gradient[0].E(0);
        double momentSq = dipole * dipole;
        if (momentSq != 0.0) {
            double s2 = 1 / r2;
            double s1 = Math.sqrt(s2);
            // normalize dr, the vector between the molecules
            drunit.E(dr);
            drunit.TE(s1);

            // v1 is the orientation of molecule 1
            IVector v1 = atom1.getOrientation().getDirection();

            // v2 is the orientation of molecule 2
            IVector v2 = atom2.getOrientation().getDirection();

            double fac = momentSq * (s2 * s1);
            double dfac = momentSq * (3 * s2 * s2 * s1);
            double udd = v1.dot(v2) - 3.0 * v1.dot(drunit) * v2.dot(drunit);
            work.Ea1Tv1(v2.dot(drunit), v1);
            work.PEa1Tv1(v1.dot(drunit), v2);
            work.PEa1Tv1(-2 * v1.dot(drunit) * v2.dot(drunit) * s1, dr);
            work.TE(3.0 * s1 * fac);
            gradient[0].PE(work);
            gradient[0].PEa1Tv1(dfac * udd, dr);

            work.E(v1);
            work.XE(v2);
            torque[0].E(v1);
            torque[0].XE(drunit);
            torque[0].TE(3.0 * v2.dot(drunit));
            torque[0].ME(work);
            torque[0].TE(fac);

            torque[1].E(v2);
            torque[1].XE(drunit);
            torque[1].TE(3.0 * v1.dot(drunit));
            torque[1].PE(work);
            torque[1].TE(fac);

            //         System.out.println("dr = " + dr);
            //         System.out.println("vi = " + v1);
            //         System.out.println("vj = " + v2);
            //         System.out.println("torque[0] = " + torque[0]);   //TODO debug only
            //         System.out.println("torque[1] = " + torque[1]);   //TODO debug only
        }

        // pairwise additive, so
        gradient[1].Ea1Tv1(-1, gradient[0]);

        return gradientAndTorque;
    }

    public Tensor[] secondDerivative(IMoleculeList molecules) {
        IMolecule molecule0 = molecules.getMolecule(0);
        IMolecule molecule1 = molecules.getMolecule(1);
        IAtomList atomList0 = molecule0.getChildList();
        IAtomList atomList1 = molecule1.getChildList();
        IAtomOriented atom0 = (IAtomOriented) atomList0.getAtom(0);
        IAtomOriented atom1 = (IAtomOriented) atomList1.getAtom(0);
        IVectorMutable pos0 = atom1.getPosition();
        IVectorMutable pos1 = atom0.getPosition();
        IVector ei = atom0.getOrientation().getDirection();
        IVector ej = atom1.getOrientation().getDirection();

        double exi = ei.getX(0);//ei and ej is the dipole orientation
        double eyi = ei.getX(1);
        double ezi = ei.getX(2);
        double exj = ej.getX(0);
        double eyj = ej.getX(1);
        double ezj = ej.getX(2);

        IVectorMutable deidxi = space.makeVector();
        IVectorMutable deidyi = space.makeVector();
        IVectorMutable deidzi = space.makeVector();
        IVectorMutable dejdxj = space.makeVector();
        IVectorMutable dejdyj = space.makeVector();
        IVectorMutable dejdzj = space.makeVector();

        double[] dejdxjD = { 0, -ezj, eyj };
        double[] dejdyjD = { ezj, 0, -exj };
        double[] dejdzjD = { -eyj, exj, 0 };
        dejdxj.E(dejdxjD);
        dejdyj.E(dejdyjD);
        dejdzj.E(dejdzjD);

        dr.Ev1Mv2(pos1, pos0);
        boundary.nearestImage(dr);//rij

        //      System.out.println("rij_P2 = " + dr);// debug only

        double r2 = dr.squared();
        secondDerivative[0].E(0);
        secondDerivative[1].E(0);
        secondDerivative[2].E(0);
        if (r2 > rCut * rCut)
            return secondDerivative;
        double r = Math.sqrt(r2);
        runit.Ea1Tv1(1.0 / r, dr);
        double coeff = dipole * dipole / r2 / r;

        //ei cross dejdxj - 3.0*(dejdxj.runit)*(ei cross runite)
        dr.E(ei);
        dr.XE(dejdxj);
        a[0].E(ei);
        a[0].XE(runit);
        a[0].TE(-3.0 * dejdxj.dot(runit));
        a[0].PE(dr);

        dr.E(ei);
        dr.XE(dejdyj);
        a[1].E(ei);
        a[1].XE(runit);
        a[1].TE(-3.0 * dejdyj.dot(runit));
        a[1].PE(dr);

        dr.E(ei);
        dr.XE(dejdzj);
        a[2].E(ei);
        a[2].XE(runit);
        a[2].TE(-3.0 * dejdzj.dot(runit));
        a[2].PE(dr);

        secondDerivative[0].E(a);
        secondDerivative[0].TE(coeff);//ij

        //      System.out.println("ei = " + ei);// debug only
        //      System.out.println("ej = " + ej);// debug only
        //      System.out.println("dudij = \n" + secondDerivative[0]);// debug only

        double[] deidxiD = { 0, -ezi, eyi };
        double[] deidyiD = { ezi, 0, -exi };
        double[] deidziD = { -eyi, exi, 0 };
        deidxi.E(deidxiD);
        deidyi.E(deidyiD);
        deidzi.E(deidziD);

        //      deidxi cross ej - 3.0*(deidxi cross runite)*(ej.runit)
        dr.E(deidxi);
        dr.XE(ej);
        a[0].E(deidxi);
        a[0].XE(runit);
        a[0].TE(-3.0 * ej.dot(runit));
        a[0].PE(dr);

        dr.E(deidyi);
        dr.XE(ej);
        a[1].E(deidyi);
        a[1].XE(runit);
        a[1].TE(-3.0 * ej.dot(runit));
        a[1].PE(dr);

        dr.E(deidzi);
        dr.XE(ej);
        a[2].E(deidzi);
        a[2].XE(runit);
        a[2].TE(-3.0 * ej.dot(runit));
        a[2].PE(dr);

        secondDerivative[1].E(a);
        secondDerivative[1].TE(coeff);//ii

        //      System.out.println("dudii = \n" + secondDerivative[1]);//TODO debug only

        //dejdxj cross ei -3.0*(ei.runit)*(dejdxj cross runit)    
        dr.E(dejdxj);
        dr.XE(ei);
        a[0].E(dejdxj);
        a[0].XE(runit);
        a[0].TE(-3.0 * ei.dot(runit));
        a[0].PE(dr);

        dr.E(dejdyj);
        dr.XE(ei);
        a[1].E(dejdyj);
        a[1].XE(runit);
        a[1].TE(-3.0 * ei.dot(runit));
        a[1].PE(dr);

        dr.E(dejdzj);
        dr.XE(ei);
        a[2].E(dejdzj);
        a[2].XE(runit);
        a[2].TE(-3.0 * ei.dot(runit));
        a[2].PE(dr);

        secondDerivative[2].E(a);
        secondDerivative[2].TE(coeff);//jj

        //      System.out.println("dudjj = \n" + secondDerivative[2]);//TODO debug only
        //      System.exit(2);// TODO bebug only
        return secondDerivative;

    }

    private static final long serialVersionUID = 1L;
    private double sigma, sigma2;
    private IBoundary boundary;
    private final IVectorMutable dr, drunit, work, runit;
    private double dipole;
    private double rCut;
    private final IVectorMutable[] gradient;
    protected final IVectorMutable[] torque;
    protected final IVectorMutable[][] gradientAndTorque;
    protected final IVectorMutable[] a;
    protected final Tensor[] secondDerivative;

}