/*
* hirshfeld.c
* newsurface
*
* Created by Anthony S. Mitchell and Joshua J. McKinnon.
* Contact: Joshua.McKinnon@une.edu.au
* Copyright (c) 2004 University of New England. All rights reserved.
*
* This file contains (almost) all of the routines required to calculate
* w(r) for the Hirshfeld surface.
*
This program is free software; you can redistribute it and/or
modify it 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.
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 General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
The GNU GPL can also be found at http://www.gnu.org
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <strings.h>
#include <unistd.h>
#include <signal.h>
#include <time.h>
#include "gdis.h"
#include "coords.h"
#include "matrix.h"
#include "space.h"
#include "numeric.h"
#include "hirshfeld.h"
#include "hirshfeld_data.h"
#include "molsurf.h"
#define MAXATOMS 54
#define NPOINTS 700
/* spline limits (in Angstroms) */
#define MAXR 7.4084815
/* crystal environment cutoff */
//#define HFS_RMAX 10.0
/* JJM DEBUG try a smaller cutoff */
#define HFS_RMAX 5.0
gboolean hfs_rho_init = FALSE;
struct rho_pak
{
gdouble x[NPOINTS];
gdouble y[NPOINTS];
gdouble y2[NPOINTS];
} atomdata[MAXATOMS];
/************************************************************/
/* initialization routine for Hirshfeld surface calculation */
/************************************************************/
void hfs_init(void)
{
gint i, j;
gdouble x /* , delta, big, small */ ;
/* only compute spline once */
if (hfs_rho_init)
return;
/* setup data array for spline fit */
for (i=0 ; i<MAXATOMS ; i++)
{
for (j=0 ; j<NPOINTS ; j++)
{
x = (gdouble) j * MAXR / (gdouble) NPOINTS;
/*** JJM: evaldens at sqrt(r) ***/
atomdata[i].x[j] = x*x;
atomdata[i].y[j] = rho[i][j];
}
/*** use a crude numerical derivative at the first point ***/
/*
delta = 0.0001;
big = evaldens(atom[i],0.0) - evaldens(atom[i],sqrt(delta));
big = big/delta;
*/
/*** use a crude numerical derivative at the last point ***/
/*
delta = 0.0001;
small = evaldens(atom[i],sqrt(maxr)) - evaldens(atom[i],(sqrt(maxr)+delta));
small = small/delta;
*/
/* FIXME - I've just used 0.0 for the end point derivatives */
spline(atomdata[i].x, atomdata[i].y, NPOINTS, 0.0, 0.0, atomdata[i].y2);
}
hfs_rho_init = TRUE;
}
/*******************************************************************/
/* return a list of cores (ie molecule) closest to the given point */
/*******************************************************************/
GSList *hfs_closest_molecule(gdouble *x, struct model_pak *model)
{
gdouble dr, min, r[3];
GSList *list;
struct core_pak *core;
struct mol_pak *mol=NULL;
/*
printf("closest to");
P3VEC(" : ", x);
*/
g_assert(model != NULL);
min = 1000000000.0;
/*
printf("min = %f\n", min);
*/
for (list=model->cores ; list ; list=g_slist_next(list))
{
core = list->data;
ARR3SET(r, core->x);
vecmat(model->latmat, r);
/*
P3VEC("core: ", r);
*/
ARR3SUB(r, x);
dr = VEC3MAGSQ(r);
if (dr < min)
{
mol = core->mol;
min = dr;
}
}
/*
printf("min = %f\n", min);
*/
if (mol)
return(mol->cores);
return(NULL);
}
/***********************************/
/* generate atoms around selection */
/***********************************/
GSList *hfs_calc_env(GSList *select, struct model_pak *model)
{
gboolean flag, primary;
gint i, imax[3];
gdouble r2, dx[3];
GSList *ilist, *list, *list2, *nlist;
struct image_pak *image;
struct core_pak *core, *copy, *test;
/* checks */
g_assert(model != NULL);
/* construct periodic image vectors */
for (i=0 ; i<model->periodic ; i++)
imax[i] = HFS_RMAX/model->pbc[i] + 1;
model->image_limit[0] = imax[0];
model->image_limit[1] = imax[0]+1;
model->image_limit[2] = imax[1];
model->image_limit[3] = imax[1]+1;
model->image_limit[4] = imax[2];
model->image_limit[5] = imax[2]+1;
space_make_images(CREATE, model);
coords_compute(model);
/* original + image iteration */
nlist = NULL;
for (list=model->cores ; list ; list=g_slist_next(list))
{
core = list->data;
/* mark cores in the central molecule */
if (g_slist_find(select, core))
primary = TRUE;
else
primary = FALSE;
/* add the atom */
ilist = NULL;
do
{
/* add coordinates */
/* NB: probably don't need to worry about shells, as this is */
/* a temporary model and the algorithm will ignore them anyway */
copy = dup_core(core);
if (ilist)
{
/* image translation */
image = ilist->data;
ARR3ADD(copy->rx, image->rx);
ilist = g_slist_next(ilist);
}
else
ilist = model->images;
/* set cartesian coords */
ARR3SET(copy->x, copy->rx);
/* JJM added line below so that copy->x contains fractional coords.
Should I take account of the centroid here? */
vecmat(model->ilatmat,copy->x);
ARR3ADD(copy->x, model->centroid);
copy->primary = primary;
/* test if core satisfies distance condition */
flag = FALSE;
for (list2=select ; list2 ; list2=g_slist_next(list2))
{
test = list2->data;
ARR3SET(dx, test->rx);
ARR3SUB(dx, copy->rx);
r2 = VEC3MAGSQ(dx);
if (r2 < HFS_RMAX*HFS_RMAX)
{
flag = TRUE;
break;
}
}
if (flag)
nlist = g_slist_prepend(nlist, copy);
}
while (ilist);
}
/* reset source model's images */
space_make_images(INITIAL, model);
#if DEBUG_HFS_CALC_ENV
printf("hfs_calc_env: %d cores\n", g_slist_length(nlist));
#endif
return(nlist);
}
/************************************/
/* retrieve an interpolated density */
/************************************/
gdouble hfs_calc_density(gint n, gdouble r2)
{
gdouble value;
/* FIXME - no support for atoms up to (and not including) Ba */
/*
g_assert(n < MAXATOMS);
*/
n = CLAMP(n, 1, MAXATOMS);
if (r2 > atomdata[n].x[NPOINTS-1])
return(0.0);
splint(atomdata[n].x, atomdata[n].y, atomdata[n].y2, NPOINTS, r2, &value);
return(value);
}
/********************************************/
/* Use the selection to calculate */
/* sensible (but not foolproof!!) limits */
/* for the marching cubes search, of the */
/* extremeties of the selection atoms in */
/* x,y,z plus <tolerance> in each direction */
/********************************************/
int hfs_calulation_limits(GSList *selection, struct model_pak *model, gdouble *min, gdouble *max)
{
gdouble tolerance[3] = {3.0, 3.0, 3.0}; /* Angstrom */
GSList *list;
struct core_pak *thisAtom;
int i;
for (i=0; i<3; i++){
min[i] = 3e10;
max[i] =-3e10;
}
vecmat(model->ilatmat,tolerance);
for (list = selection; list; list=g_slist_next(list)){
thisAtom = list->data;
for (i=0; i<3; i++){
if ( thisAtom->x[i] < min[i]) min[i] = thisAtom->x[i];
if ( thisAtom->x[i] > max[i]) max[i] = thisAtom->x[i];
}
}
for (i=0; i<3; i++) {
min[i] = min[i] - tolerance[i];
max[i] = max[i] + tolerance[i];
}
return(0);
}
/*************************************/
/* Use the gradient of the Hirshfeld */
/* surface to calculate the normal */
/*************************************/
int hfs_calc_normals(GSList *points, struct model_pak *model, GSList *myThing)
{
gdouble delta = 5.0e-3;
gdouble gx, gy, gz;
gdouble wPoint, wgx, wgy, wgz;
gdouble xDelta[3],yDelta[3],zDelta[3];
gdouble temp[3];
GSList *pointList;
struct smv_pak *thisPoint;
time_t jobStartTime, jobEndTime;
jobStartTime = time(NULL);
delta = 5.0e-3;
for (pointList=points; pointList; pointList=g_slist_next(pointList)){
thisPoint = pointList->data;
if (!thisPoint->adj)
continue;
ARR3SET(temp, thisPoint->rx);
gx = thisPoint->rx[0] + delta;
gy = thisPoint->rx[1] + delta;
gz = thisPoint->rx[2] + delta;
VEC3SET(xDelta,gx,thisPoint->rx[1],thisPoint->rx[2]);
VEC3SET(yDelta,thisPoint->rx[0],gy,thisPoint->rx[2]);
VEC3SET(zDelta,thisPoint->rx[0],thisPoint->rx[1],gz);
wPoint = hfs_calc_wf(thisPoint->rx[0],thisPoint->rx[1],thisPoint->rx[2],model,myThing);
wgx = hfs_calc_wf(xDelta[0],xDelta[1],xDelta[2],model,myThing);
wgy = hfs_calc_wf(yDelta[0],yDelta[1],yDelta[2],model,myThing);
wgz = hfs_calc_wf(zDelta[0],zDelta[1],zDelta[2],model,myThing);
thisPoint->nx[0] = (wgx - wPoint)/delta;
thisPoint->nx[1] = (wgy - wPoint)/delta;
thisPoint->nx[2] = (wgz - wPoint)/delta;
// normalize(thisPoint->nx,3);
ARR3SET(thisPoint->n, thisPoint->nx);
// P3VEC("Normal:",thisPoint->nx);
}
jobEndTime = time(NULL);
#if DEBUG_HFS_NORMALS
fprintf(stderr," [%.1fs]\n",(float)(jobEndTime - jobStartTime));
#endif
return(0);
}
/*************************************/
/* Use the gradient of the promolecule */
/* surface to calculate the normal */
/*************************************/
int ssatoms_calc_normals(GSList *points, struct model_pak *model)
{
gdouble delta = 5.0e-3;
gdouble gx, gy, gz;
gdouble wPoint, wgx, wgy, wgz;
gdouble xDelta[3],yDelta[3],zDelta[3];
gdouble temp[3];
GSList *pointList;
struct smv_pak *thisPoint;
time_t jobStartTime, jobEndTime;
jobStartTime = time(NULL);
delta = 5.0e-3;
for (pointList=points; pointList; pointList=g_slist_next(pointList)){
thisPoint = pointList->data;
if (!thisPoint->adj)
continue;
ARR3SET(temp, thisPoint->rx);
gx = thisPoint->rx[0] + delta;
gy = thisPoint->rx[1] + delta;
gz = thisPoint->rx[2] + delta;
VEC3SET(xDelta,gx,thisPoint->rx[1],thisPoint->rx[2]);
VEC3SET(yDelta,thisPoint->rx[0],gy,thisPoint->rx[2]);
VEC3SET(zDelta,thisPoint->rx[0],thisPoint->rx[1],gz);
wPoint = ssatoms_calc_density(thisPoint->rx[0],thisPoint->rx[1],thisPoint->rx[2],model);
wgx = ssatoms_calc_density(xDelta[0],xDelta[1],xDelta[2],model);
wgy = ssatoms_calc_density(yDelta[0],yDelta[1],yDelta[2],model);
wgz = ssatoms_calc_density(zDelta[0],zDelta[1],zDelta[2],model);
thisPoint->nx[0] = (wgx - wPoint)/delta;
thisPoint->nx[1] = (wgy - wPoint)/delta;
thisPoint->nx[2] = (wgz - wPoint)/delta;
normalize(thisPoint->nx,3);
ARR3SET(thisPoint->n, thisPoint->nx);
// P3VEC("Normal:",thisPoint->nx);
}
jobEndTime = time(NULL);
#if DEBUG_HFS_NORMALS
fprintf(stderr," [%.1fs]\n",(float)(jobEndTime - jobStartTime));
#endif
return(0);
}
/******************************************/
/* compute the weight function at a point */
/* */
/* Note we now pass three doubles (x,y,z) */
/* instead of a single vector because it */
/* makes life much easier for the second */
/* derivative calculations */
/******************************************/
gdouble hfs_calc_wf(gdouble x, gdouble y, gdouble z, struct model_pak *model, GSList *myThing)
{
gdouble r2;
gdouble weight_r=0.0, r[3];
gdouble point[3]; /*vector for the point (x,y,z) */
gdouble rhoNugget, rhoMol, rhoXtal;
GSList *list, *mlist;
struct core_pak *core;
VEC3SET(point,x,y,z);
/* checks */
g_assert(model != NULL);
g_assert(hfs_rho_init == TRUE);
/* add up the density contribution from the molecule */
/*** JJM The 'molecule' that we're calculating the Hirshfeld surface of should be the selected molecule. Hence the new definition of mlist here ***/
mlist = model->selection;
rhoMol=0.0;
/***J loop over the atoms in the selected molecule J***/
for (list=mlist ; list ; list=g_slist_next(list))
{
core = list->data;
/* compute cartesian distance squared (in Angs) */
// vecmat(model->latmat,core->x);
ARR3SET(r, core->rx);
vecmat(model->ilatmat,r);
ARR3ADD(r, model->centroid);
vecmat(model->latmat,r);
ARR3SUB(r, point);
r2 = VEC3MAGSQ(r);
/* ARR3SET(coreCartesian, core->x);
vecmat(model->latmat,coreCartesian);
ARR3SUB(coreCartesian,x);
r2 = VEC3MAGSQ(r);
*/
rhoNugget = hfs_calc_density(core->atom_code, r2);
rhoMol += rhoNugget;
}
/* add up the total density */
rhoXtal=0.0;
/***JJM loop over the atoms in the crystal J***/
for (list=myThing ; list ; list=g_slist_next(list))
{
core = list->data;
/* compute cartesian distance squared (in Angs) */
// vecmat(model->latmat,core->x);
ARR3SET(r, core->rx);
vecmat(model->ilatmat,r);
ARR3ADD(r, model->centroid);
vecmat(model->latmat,r);
#ifdef JJM_DEBUG_EXTREME02
P3VEC(" r: ", r);
P3VEC(" x: ", point);
#endif
ARR3SUB(r, point);
r2 = VEC3MAGSQ(r);
/* ARR3SET(coreCartesian, core->x);
vecmat(model->latmat,coreCartesian);
ARR3SUB(coreCartesian,x);
r2 = VEC3MAGSQ(r);
*/
rhoNugget = hfs_calc_density(core->atom_code, r2);
rhoXtal += rhoNugget;
}
#ifdef JJM_DEBUG02
if ( rhoMol > 0.0 || rhoXtal > 0.0){
fprintf(stderr,"rhoMol: %.3f ___ rhoXtal: %.3f\n",rhoMol,rhoXtal);
}
#endif
/* calc ratio */
if (rhoXtal != 0.0){
/* NEW */
if (rhoXtal > rhoMol)
weight_r = rhoMol/rhoXtal;
else {
fprintf(stderr,"WARNING: Weight function > 1.0. Value: %f\n",weight_r);
weight_r = 0.0;
}
}
else
weight_r = 0.0;
#define STRICT 1
#if STRICT
g_assert(weight_r <= 1.0);
#else
if (weight_r > 1.0)
weight_r = 1.0;
#endif
/* JJM debug
if (weight_r > 0.5){
fprintf(stderr,"w %.2f ",weight_r);
P3VEC("at ",x);
}
*/
return(weight_r);
}
/*********************************************************************/
/* compute the sum of spherical atom density at a point */
/* a.k.a. Promolecule */
/* See Mitchell & Spackman, J. Comp. Chem., v. 21, pp 933-942 (2000) */
/* */
/* */
/* */
/* Note we now pass three doubles (x,y,z) */
/* instead of a single vector because it */
/* makes life much easier for the second */
/* derivative calculations */
/*********************************************************************/
gdouble ssatoms_calc_density(gdouble x, gdouble y, gdouble z, struct model_pak *model)
{
gdouble r2;
gdouble r[3];
gdouble point[3]; /*vector for the point (x,y,z) */
gdouble rhoNugget, rhoMol;
GSList *list, *mlist;
struct core_pak *core;
VEC3SET(point,x,y,z);
/* checks */
g_assert(model != NULL);
g_assert(hfs_rho_init == TRUE);
/* add up the density contribution from the molecule */
/*** JJM The 'molecule' that we're calculating the surface of should be the selected molecule. Hence the new definition of mlist here ***/
mlist = model->selection;
rhoMol=0.0;
/***J loop over the atoms in the selected molecule J***/
for (list=mlist ; list ; list=g_slist_next(list))
{
core = list->data;
/* compute cartesian distance squared (in Angs) */
// vecmat(model->latmat,core->x);
ARR3SET(r, core->rx);
vecmat(model->ilatmat,r);
ARR3ADD(r, model->centroid);
vecmat(model->latmat,r);
ARR3SUB(r, point);
r2 = VEC3MAGSQ(r);
/* ARR3SET(coreCartesian, core->x);
vecmat(model->latmat,coreCartesian);
ARR3SUB(coreCartesian,x);
r2 = VEC3MAGSQ(r);
*/
rhoNugget = hfs_calc_density(core->atom_code, r2);
rhoMol += rhoNugget;
}
#ifdef JJM_DEBUG02
if ( rhoMol > 0.0 ){
fprintf(stderr,"rhoMol: %.3f\n",rhoMol);
}
#endif
return(rhoMol);
}
syntax highlighted by Code2HTML, v. 0.9.1