Package stallone.potential

Source Code of stallone.potential.CoulombSystem

/*
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
package stallone.potential;

import stallone.api.doubles.IDifferentiableMetric;
import stallone.api.doubles.Doubles;
import stallone.api.doubles.IDoubleArray;
import stallone.api.ints.IIntArray;


/**
*
* @author noe
*/
public class CoulombSystem  extends AbstractPotential
{
    private static double epsilon0 = 0.0005727657013; // in elementary charges * nm^-1 * kJ mol^-1
    private double epsilonr = 1;

    private IDoubleArray coordinates;
    private IDifferentiableMetric<IDoubleArray> metric;
    private int natoms = 0;
    private boolean[][] bonded;
    // particle radius
    private IDoubleArray charges;

    private double energy;
    private IDoubleArray gradient;

    public CoulombSystem(IDoubleArray _charges, double _epsilonr)
    {
        this.natoms = _charges.size();
        this.bonded = new boolean[_charges.size()][_charges.size()];
        this.charges = _charges;
        this.epsilonr = _epsilonr;

        this.gradient = Doubles.create.array(_charges.size(), 3);
    }

    public CoulombSystem(IDoubleArray _charges, IIntArray _bonded, double _epsilonr)
    {
        if (_bonded.columns() != 2)
            throw(new IllegalArgumentException("_bonded must be a nx2 table"));

        this.natoms = _charges.size();
        this.bonded = new boolean[_charges.size()][_charges.size()];
        for (int i=0; i<_bonded.size(); i++)
            this.bonded[_bonded.get(i,0)][_bonded.get(i,1)] = true;
        this.charges = _charges;
        this.epsilonr = _epsilonr;

        this.gradient = Doubles.create.array(_charges.size(), 3);
    }

    public CoulombSystem(IDoubleArray _charges, boolean[][] _bonded, double _epsilonr)
    {
        this.natoms = _charges.size();
        this.bonded = _bonded;
        this.charges = _charges;
        this.epsilonr = _epsilonr;

        this.gradient = Doubles.create.array(_charges.size(), 3);
    }

    public void setMetric(IDifferentiableMetric<IDoubleArray> m)
    {
        metric = m;
    }

    @Override
    public void setCoordinates(IDoubleArray _coordinates)
    {
        coordinates = (IDoubleArray) _coordinates.copy();
    }

    @Override
    public boolean calculate()
    {
        // initialize energy and gradient
        this.energy = 0;
        for (int i = 0; i < gradient.size(); i++)
        {
            this.gradient.set(i, 0);
        }

        // go through all pairs
        for (int i = 0; i < natoms - 1; i++)
        {
            for (int j = i + 1; j < natoms; j++)
            {
                if (bonded[i][j])
                {
                    continue;
                }

                // get actual distance
                IDoubleArray c1 = coordinates.viewRow(i);
                IDoubleArray c2 = coordinates.viewRow(j);
                double d = metric.distance(c1, c2);

                // energy
                double e = charges.get(i) * charges.get(j) / (4*Math.PI*epsilon0*epsilonr*d);
                this.energy += e;
//System.out.println(" c ["+i+","+j+"]\t"+e);

                // gradient
                double fm = -charges.get(i) * charges.get(j) / (4*Math.PI*epsilon0*epsilonr*d*d);

                IDoubleArray mg1 = metric.gradientX(c1, c2);
                for (int dim = 0; dim < 3; dim++)
                {
                    this.gradient.set(i, dim, this.gradient.get(i, dim) + fm * mg1.get(dim));
                }

                IDoubleArray mg2 = metric.gradientY(c1, c2);
                for (int dim = 0; dim < 3; dim++)
                {
                    this.gradient.set(j, dim, this.gradient.get(j, dim) + fm * mg2.get(dim));
                }
            }
        }

        return (true);
    }

    @Override
    public int getNDimensions()
    {
        return (natoms * 3);
    }

    @Override
    public double getEnergy()
    {
        return (energy);
    }

    @Override
    public IDoubleArray getGradient()
    {
        return (gradient);
    }

    @Override
    public IDoubleArray getCoordinates()
    {
        return (coordinates);
    }

        @Override
    public int getNumberOfVariables()
    {
        return(coordinates.size());
    }


}
TOP

Related Classes of stallone.potential.CoulombSystem

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.