Commit f3a81bc3 by Francois Gygi

Merge branch 'develop' into vext

parents 6fe308c2 a23e6450
This diff is collapsed. Click to expand it.
Qbox First-Principles Molecular Dynamics
========================================
Qbox is a C++/MPI/OpenMP implementation of First-Principles Molecular Dynamics.
It implements electronic structure calculations within the framework
of Density Functional Theory, using a plane-wave basis set and
norm-conserving pseudopotentials.
For documentation, examples and installation instructions
see http://qboxcode.org
Qbox is developed in the Gygi group at the University of California Davis.
Qbox is distributed under the terms of the GNU General Public License
as published by the Free Software Foundation, either version 2 of
the License, or (at your option) any later version.
See the file COPYING in the root directory of this distribution
or <http://www.gnu.org/licenses/>.
Qbox was written by Francois Gygi.
Ivan Duchemin developed the Hartree-Fock exchange implementation and the
client-server interface.
Jun Wu developed the implementation of van der Waals density functionals.
Quan (Andy) Wan implemented the FFTW3 interface.
William Dawson contributed to the implementation of Hartree-Fock exchange.
Martin Schlipf implemented Optimized Norm-Conserving Vanderbilt potentials.
Martin Schlipf implemented the HSE hybrid density functional
Quan (Andy) Wan implemented the electric field.
Ma He contributed to the implementation of the response functionality.
#-------------------------------------------------------------------------------
#
# Copyright (c) 2008-2017 The Regents of the University of California
#
# This file is part of Qbox
#
# Qbox is distributed under the terms of the GNU General Public License
# as published by the Free Software Foundation, either version 2 of
# the License, or (at your option) any later version.
# See the file COPYING in the root directory of this distribution
# or <http://www.gnu.org/licenses/>.
#
#-------------------------------------------------------------------------------
#
# centos7.mk
#
#-------------------------------------------------------------------------------
# Prerequisites:
# On a Centos-7.3 system, install the following packages:
# yum install xerces-c xerces-c-devel
# yum install openmpi openmpi-devel
# yum install lapack lapack-devel
# yum install fftw fftw-devel
# yum install scalapack-common scalapack-openmpi \
# scalapack-openmpi-devel scalapack-openmpi-static
# yum install libuuid libuuid-devel
#
#-------------------------------------------------------------------------------
PLT=Linux_x8664
#-------------------------------------------------------------------------------
PLTOBJECTS = readTSC.o
CXX=mpicxx
LD=$(CXX)
PLTFLAGS += -DIA32 -D_LARGEFILE_SOURCE \
-D_FILE_OFFSET_BITS=64 -DUSE_MPI -DSCALAPACK -DADD_ \
-DAPP_NO_THREADS -DXML_USE_NO_THREADS -DUSE_XERCES \
-DXERCESC_3 -DMPICH_IGNORE_CXX_SEEK -DUSE_UUID
# FFT must be FFTW2, FFTW3, ESSL or NOLIB
FFT=FFTW3
ifeq ($(FFT),FFTW2)
PLTFLAGS += -DUSE_FFTW2
#PLTFLAGS += -DFFTWMEASURE
#FFTWDIR=$(HOME)/software/fftw/Linux_x8664/fftw-2.1.5/fftw
#FFTWINCLUDEDIR=$(FFTWDIR)
#FFTWLIBDIR=$(FFTWDIR)/.libs
#INCLUDE += -I$(FFTWINCLUDEDIR)
#LIBPATH += -L$(FFTWLIBDIR)
LIBS += -lfftw
endif
ifeq ($(FFT),FFTW3)
PLTFLAGS += -DUSE_FFTW3
#PLTFLAGS += -DFFTWMEASURE
#PLTFLAGS += -DFFTW_TRANSPOSE
PLTFLAGS += -DFFTW3_2D
#FFTWDIR=$(HOME)/software/fftw/fftw-3.3.4
#FFTWINCLUDEDIR=$(FFTWDIR)/api
#FFTWLIBDIR=$(FFTWDIR)/.libs
#INCLUDE += -I$(FFTWINCLUDEDIR)
#LIBPATH += -L$(FFTWLIBDIR)
LIBS += -lfftw3
endif
ifeq ($(FFT),ESSL)
$(error ESSL library not available)
endif
ifeq ($(FFT),NOLIB)
PLTFLAGS += -DFFT_NOLIB
endif
CXXFLAGS= -g -O3 -Wunused -D$(PLT) $(INCLUDE) $(PLTFLAGS) $(DFLAGS)
LIBS += -lpthread -lxerces-c -lscalapack -llapack -lblas -luuid
LDFLAGS = $(LIBPATH) $(LIBS)
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#
# Copyright (c) 2008-2018 The Regents of the University of California
#
# This file is part of Qbox
#
# Qbox is distributed under the terms of the GNU General Public License
# as published by the Free Software Foundation, either version 2 of
# the License, or (at your option) any later version.
# See the file COPYING in the root directory of this distribution
# or <http://www.gnu.org/licenses/>.
#
#-------------------------------------------------------------------------------
#
# centos7_intel.mk
#
#-------------------------------------------------------------------------------
# build on CentOS 7 using the intel compiler and the threaded MKL library
# with openmpi
#
# Prerequisites:
# On a Centos-7.3 system, install the following packages:
# yum install xerces-c xerces-c-devel
# yum install openmpi openmpi-devel
# yum install fftw fftw-devel
# yum install libuuid libuuid-devel
# intel compiler and MKL library
#
#-------------------------------------------------------------------------------
# Use the following definitions to compile and use qb
# export LD_LIBRARY_PATH=/opt/intel/mkl/intel64/lib
# export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/intel/lib/intel64
# export PATH=$PATH:/opt/intel/bin
MKLROOT = /opt/intel/mkl
PLTOBJECTS = readTSC.o
CXX=env OMPI_CXX=icc mpicxx
LD=$(CXX)
PLTFLAGS += -D_LARGEFILE_SOURCE \
-D_FILE_OFFSET_BITS=64 -DUSE_MPI -DSCALAPACK -DADD_ \
-DAPP_NO_THREADS -DXML_USE_NO_THREADS -DUSE_XERCES \
-DXERCESC_3 -DMPICH_IGNORE_CXX_SEEK -DUSE_UUID \
-DUSE_FFTW3 -DFFTW3_2D
LIBS = -lfftw3
INCLUDE = -I$(MKLROOT)/include
CXXFLAGS= -g -O3 -openmp $(INCLUDE) $(PLTFLAGS)
LIBPATH = -L$(MKLROOT)/lib/intel64
LIBS += -Wl,-Bstatic -Wl,--start-group \
-lmkl_scalapack_lp64 -lmkl_blacs_openmpi_lp64 \
-lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread \
-Wl,--end-group -Wl,-Bdynamic \
-lxerces-c -luuid -liomp5
LDFLAGS = $(LIBPATH) $(LIBS)
#-------------------------------------------------------------------------------
#
# Copyright (c) 2014 The Regents of the University of California
#
# This file is part of Qbox
#
# Qbox is distributed under the terms of the GNU General Public License
# as published by the Free Software Foundation, either version 2 of
# the License, or (at your option) any later version.
# See the file COPYING in the root directory of this distribution
# or <http://www.gnu.org/licenses/>.
#
#-------------------------------------------------------------------------------
#
# macosx_openmpi.mk
#
#-------------------------------------------------------------------------------
#
PLT=MacOSX_x86_64
#-------------------------------------------------------------------------------
SCALAPACKDIR=$(HOME)/software/scalapack/scalapack-2.0.2
LAPACKDIR=$(HOME)/software/lapack/lapack-3.5.0
BLASDIR=$(HOME)/software/blas/BLAS
PLTOBJECTS = readTSC.o
CXX=mpiCC
LD=$(CXX)
PLTFLAGS += -DIA32 -DUSE_MPI -DUSE_FFTW2 -D_LARGEFILE_SOURCE \
-D_FILE_OFFSET_BITS=64 -DADD_ \
-DAPP_NO_THREADS -DXML_USE_NO_THREADS -DUSE_XERCES -DXERCESC_3 \
-DSCALAPACK -DUSE_UUID
INCLUDE = -I$(FFTWDIR)
CXXFLAGS= -g -O3 -Wunused -D$(PLT) $(INCLUDE) $(PLTFLAGS) $(DFLAGS)
LIBPATH = -L$(SCALAPACKDIR) -L$(LAPACKDIR) -L$(BLASDIR)
LIBS = -lfftw -lscalapack -llapack -lblas -lm -lgfortran \
-lxerces-c -lpthread -lmpi_cxx -lstdc++ -lcurl -liconv
LDFLAGS = $(LIBPATH) $(LIBS)
#-------------------------------------------------------------------------------
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2014 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// AlphaRSH.h
//
////////////////////////////////////////////////////////////////////////////////
#ifndef ALPHARSH_H
#define ALPHARSH_H
#include<iostream>
#include<iomanip>
#include<sstream>
#include<stdlib.h>
#include "Sample.h"
class AlphaRSH : public Var
{
Sample *s;
public:
const char *name ( void ) const { return "alpha_RSH"; };
int set ( int argc, char **argv )
{
if ( argc != 2 )
{
if ( ui->onpe0() )
cout << " alpha_RSH takes only one value" << endl;
return 1;
}
double v = atof(argv[1]);
if ( v < 0.0 )
{
if ( ui->onpe0() )
cout << " alpha_RSH must be non-negative" << endl;
return 1;
}
s->ctrl.alpha_RSH = v;
return 0;
}
string print (void) const
{
ostringstream st;
st.setf(ios::left,ios::adjustfield);
st << setw(10) << name() << " = ";
st.setf(ios::right,ios::adjustfield);
st << s->ctrl.alpha_RSH;
return st.str();
}
AlphaRSH(Sample *sample) : s(sample)
{
s->ctrl.alpha_RSH = 0;
}
};
#endif
......@@ -80,7 +80,7 @@ vector<vector<double> > &rp) const
D3vector r3(pr3);
D3vector g1,g2,g3;
grad_sigma(r1,r2,r3,g1,g2,g3);
const double a = bond_angle(r1,r2,r3);
// const double a = bond_angle(r1,r2,r3);
double ng = g1*g1 + g2*g2 + g3*g3;
assert(ng>=0.0);
......@@ -358,9 +358,9 @@ void AngleConstraint::grad_sigma(const D3vector &r1, const D3vector &r2,
// angle is large enough. Use finite differences
//cout << " ========= grad_sigma using finite differences" << endl;
const double r12_inv = 1.0/length(r12);
const double r32_inv = 1.0/length(r32);
const double a = bond_angle(r1,r2,r3);
// const double r12_inv = 1.0/length(r12);
// const double r32_inv = 1.0/length(r32);
// const double a = bond_angle(r1,r2,r3);
const double l12 = length(r1-r2);
const double l32 = length(r3-r2);
......
......@@ -503,6 +503,25 @@ void AtomSet::randomize_velocities(double temp)
}
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::rcm(void) const
{
D3vector mrsum;
double msum = 0.0;
for ( int is = 0; is < atom_list.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < atom_list[is].size(); ia++ )
{
D3vector r = atom_list[is][ia]->position();
mrsum += mass * r;
msum += mass;
}
}
if ( msum == 0.0 ) return D3vector(0,0,0);
return mrsum / msum;
}
////////////////////////////////////////////////////////////////////////////////
D3vector AtomSet::vcm(void) const
{
D3vector mvsum;
......@@ -542,6 +561,141 @@ void AtomSet::reset_vcm(void)
}
////////////////////////////////////////////////////////////////////////////////
void AtomSet::reset_rotation(void)
{
// check for special case of zero or one atom
if ( size() < 2 ) return;
D3vector rc = rcm();
D3vector vc = vcm();
vector<vector<double> > rt;
get_positions(rt);
vector<vector<double> > vt;
get_velocities(vt);
// compute angular momentum w.r.t. the center of mass
D3vector L;
for ( int is = 0; is < vt.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
D3vector v(vt[is][3*ia+0],vt[is][3*ia+1],vt[is][3*ia+2]);
v -= vc;
L += mass * ( r ^ v );
}
}
// compute inertia tensor a and vector omega
// check for special case of all atoms aligned, for which the
// inertia tensor has a zero eigenvalue
// check if all atoms are aligned
// collect all positions in a single vector of D3vectors
vector<D3vector> rv(size());
int iat = 0;
for ( int is = 0; is < vt.size(); is++ )
for ( int ia = 0; ia < na(is); ia++ )
rv[iat++] = D3vector(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
// normalized direction e = (rv[1]-rv[0])
D3vector e = normalized(rv[1]-rv[0]);
bool aligned = true;
for ( int i = 2; (i < size()) && aligned; i++ )
{
D3vector u = normalized(rv[i]-rv[0]);
aligned &= length(u^e) < 1.e-6;
}
D3vector omega;
// compute the inertia tensor
// treat separately the case of all atoms aligned
if ( aligned )
{
// inertia tensor reduces to a scalar
double a = 0.0;
for ( int is = 0; is < vt.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
a += mass * norm2(r);
}
}
omega = L / a;
}
else
{
double a00,a01,a02,a10,a11,a12,a20,a21,a22;
a00=a01=a02=a10=a11=a12=a20=a21=a22=0.0;
for ( int is = 0; is < vt.size(); is++ )
{
double mass = species_list[is]->mass();
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
a00 += mass * ( r.y * r.y + r.z * r.z );
a11 += mass * ( r.x * r.x + r.z * r.z );
a22 += mass * ( r.x * r.x + r.y * r.y );
a01 -= mass * r.x * r.y;
a02 -= mass * r.x * r.z;
a12 -= mass * r.y * r.z;
}
}
a10 = a01;
a20 = a02;
a21 = a12;
// inverse b of the inertia tensor a
// determinant
double det = a00 * ( a11 * a22 - a21 * a12 ) -
a01 * ( a10 * a22 - a20 * a12 ) +
a02 * ( a10 * a21 - a20 * a11 );
// the determinant must be positive
assert(det>1.e-8);
double b00,b01,b02,b10,b11,b12,b20,b21,b22;
b00 = ( -a12*a21 + a11*a22 ) / det;
b10 = ( a12*a20 - a10*a22 ) / det;
b20 = ( -a11*a20 + a10*a21 ) / det;
b01 = ( a02*a21 - a01*a22 ) / det;
b11 = ( -a02*a20 + a00*a22 ) / det;
b21 = ( a01*a20 - a00*a21 ) / det;
b02 = ( -a02*a11 + a01*a12 ) / det;
b12 = ( a02*a10 - a00*a12 ) / det;
b22 = ( -a01*a10 + a00*a11 ) / det;
// omega = inverse(a) * L
omega.x = b00 * L.x + b01 * L.y + b02 * L.z;
omega.y = b10 * L.x + b11 * L.y + b12 * L.z;
omega.z = b20 * L.x + b21 * L.y + b22 * L.z;
}
// correct velocities: v = v - omega ^ r
for ( int is = 0; is < vt.size(); is++ )
{
for ( int ia = 0; ia < na(is); ia++ )
{
D3vector r(rt[is][3*ia+0],rt[is][3*ia+1],rt[is][3*ia+2]);
r -= rc;
D3vector dv = omega ^ r;
vt[is][3*ia+0] -= dv.x;
vt[is][3*ia+1] -= dv.y;
vt[is][3*ia+2] -= dv.z;
}
}
set_velocities(vt);
}
////////////////////////////////////////////////////////////////////////////////
void AtomSet::fold_in_ws(void)
{
vector<vector<double> > p;
......
......@@ -85,10 +85,12 @@ class AtomSet
void rescale_velocities(double fac);
void randomize_velocities(double temp);
void randomize_positions(double amplitude);
D3vector rcm(void) const;
D3vector vcm(void) const;
D3vector dipole(void) const;
D3tensor quadrupole(void) const;
void reset_vcm(void);
void reset_rotation(void);
void fold_in_ws(void);
int size(void) const;
};
......
......@@ -42,7 +42,7 @@ void AtomSetHandler::startElement(const XMLCh* const uri,
const Attributes& attributes)
{
// cout << " AtomSetHandler::startElement " << StrX(qname) << endl;
string locname(XMLString::transcode(localname));
string locname = StrX(localname).localForm();
// consider only elements that are dealt with directly by AtomSetHandler
// i.e. "atom". The "species" element is delegated to a SpeciesHandler
......@@ -52,8 +52,8 @@ void AtomSetHandler::startElement(const XMLCh* const uri,
unsigned int len = attributes.getLength();
for (unsigned int index = 0; index < len; index++)
{
string attrname(XMLString::transcode(attributes.getLocalName(index)));
string attrval(XMLString::transcode(attributes.getValue(index)));
string attrname = StrX(attributes.getLocalName(index)).localForm();
string attrval = StrX(attributes.getValue(index)).localForm();
istringstream stst(attrval);
if ( attrname == "a")
{
......@@ -79,7 +79,7 @@ void AtomSetHandler::startElement(const XMLCh* const uri,
unsigned int len = attributes.getLength();
for (unsigned int index = 0; index < len; index++)
{
string attrname(XMLString::transcode(attributes.getLocalName(index)));
string attrname = StrX(attributes.getLocalName(index)).localForm();
if ( attrname == "name")
{
current_atom_name = StrX(attributes.getValue(index)).localForm();
......@@ -96,7 +96,7 @@ void AtomSetHandler::startElement(const XMLCh* const uri,
void AtomSetHandler::endElement(const XMLCh* const uri,
const XMLCh* const localname, const XMLCh* const qname, string& content)
{
string locname(XMLString::transcode(localname));
string locname = StrX(localname).localForm();
// cout << " AtomSetHandler::endElement " << locname << endl;
istringstream stst(content);
if ( locname == "unit_cell")
......@@ -146,14 +146,14 @@ StructureHandler* AtomSetHandler::startSubHandler(const XMLCh* const uri,
// If it can, return a pointer to the StructureHandler, otherwise return 0
// cout << " AtomSetHandler::startSubHandler " << StrX(qname) << endl;
string locname(XMLString::transcode(localname));
string locname = StrX(localname).localForm();
if ( locname == "species")
{
// check for species attributes
unsigned int len = attributes.getLength();
for (unsigned int index = 0; index < len; index++)
{
string attrname(XMLString::transcode(attributes.getLocalName(index)));
string attrname = StrX(attributes.getLocalName(index)).localForm();
if ( attrname == "name")
{
current_species_name = StrX(attributes.getValue(index)).localForm();
......@@ -175,7 +175,7 @@ void AtomSetHandler::endSubHandler(const XMLCh* const uri,
const XMLCh* const localname, const XMLCh* const qname,
const StructureHandler* const last)
{
string locname(XMLString::transcode(localname));
string locname = StrX(localname).localForm();
// cout << " AtomSetHandler::endSubHandler " << locname << endl;
if ( locname == "species" )
{
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// BHandHLYPFunctional.C
//
////////////////////////////////////////////////////////////////////////////////
#include <cmath>
#include <cassert>
#include "BHandHLYPFunctional.h"
#include "BLYPFunctional.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
BHandHLYPFunctional::BHandHLYPFunctional(const vector<vector<double> > &rhoe)
{
_nspin = rhoe.size();
if ( _nspin > 1 ) assert(rhoe[0].size() == rhoe[1].size());
_np = rhoe[0].size();
if ( _nspin == 1 )
{
_exc.resize(_np);
_vxc1.resize(_np);
_vxc2.resize(_np);
_grad_rho[0].resize(_np);
_grad_rho[1].resize(_np);
_grad_rho[2].resize(_np);
rho = &rhoe[0][0];
grad_rho[0] = &_grad_rho[0][0];
grad_rho[1] = &_grad_rho[1][0];
grad_rho[2] = &_grad_rho[2][0];
exc = &_exc[0];
vxc1 = &_vxc1[0];
vxc2 = &_vxc2[0];
}
else
{
_exc_up.resize(_np);
_exc_dn.resize(_np);
_vxc1_up.resize(_np);
_vxc1_dn.resize(_np);
_vxc2_upup.resize(_np);
_vxc2_updn.resize(_np);
_vxc2_dnup.resize(_np);
_vxc2_dndn.resize(_np);
_grad_rho_up[0].resize(_np);
_grad_rho_up[1].resize(_np);
_grad_rho_up[2].resize(_np);
_grad_rho_dn[0].resize(_np);
_grad_rho_dn[1].resize(_np);
_grad_rho_dn[2].resize(_np);
rho_up = &rhoe[0][0];
rho_dn = &rhoe[1][0];
grad_rho_up[0] = &_grad_rho_up[0][0];
grad_rho_up[1] = &_grad_rho_up[1][0];
grad_rho_up[2] = &_grad_rho_up[2][0];
grad_rho_dn[0] = &_grad_rho_dn[0][0];
grad_rho_dn[1] = &_grad_rho_dn[1][0];
grad_rho_dn[2] = &_grad_rho_dn[2][0];
exc_up = &_exc_up[0];
exc_dn = &_exc_dn[0];
vxc1_up = &_vxc1_up[0];
vxc1_dn = &_vxc1_dn[0];
vxc2_upup = &_vxc2_upup[0];
vxc2_updn = &_vxc2_updn[0];
vxc2_dnup = &_vxc2_dnup[0];
vxc2_dndn = &_vxc2_dndn[0];
}
}
////////////////////////////////////////////////////////////////////////////////
void BHandHLYPFunctional::setxc(void)
{
if ( _np == 0 ) return;
if ( _nspin == 1 )
{
assert( rho != 0 );
assert( grad_rho[0] != 0 && grad_rho[1] != 0 && grad_rho[2] != 0 );
assert( exc != 0 );
assert( vxc1 != 0 );
assert( vxc2 != 0 );
for ( int i = 0; i < _np; i++ )
{
double grad = sqrt(grad_rho[0][i]*grad_rho[0][i] +
grad_rho[1][i]*grad_rho[1][i] +
grad_rho[2][i]*grad_rho[2][i] );
excbhandhlyp(rho[i],grad,&exc[i],&vxc1[i],&vxc2[i]);
}
}
else
{
assert( rho_up != 0 );
assert( rho_dn != 0 );
assert( grad_rho_up[0] != 0 && grad_rho_up[1] != 0 && grad_rho_up[2] != 0 );
assert( grad_rho_dn[0] != 0 && grad_rho_dn[1] != 0 && grad_rho_dn[2] != 0 );
assert( exc_up != 0 );
assert( exc_dn != 0 );
assert( vxc1_up != 0 );
assert( vxc1_dn != 0 );
assert( vxc2_upup != 0 );
assert( vxc2_updn != 0 );
assert( vxc2_dnup != 0 );
assert( vxc2_dndn != 0 );
for ( int i = 0; i < _np; i++ )
{
double grx_up = grad_rho_up[0][i];
double gry_up = grad_rho_up[1][i];
double grz_up = grad_rho_up[2][i];
double grx_dn = grad_rho_dn[0][i];
double gry_dn = grad_rho_dn[1][i];
double grz_dn = grad_rho_dn[2][i];
double grad_up2 = grx_up*grx_up + gry_up*gry_up + grz_up*grz_up;
double grad_dn2 = grx_dn*grx_dn + gry_dn*gry_dn + grz_dn*grz_dn;
double grad_up_grad_dn = grx_up*grx_dn + gry_up*gry_dn + grz_up*grz_dn;