| 1 |
#define EXTERN extern |
| 2 |
|
| 3 |
#include "pcwin.h" |
| 4 |
#include "utility.h" |
| 5 |
|
| 6 |
EXTERN struct t_minim_values { |
| 7 |
int iprint, ndc, nconst; |
| 8 |
float dielc; |
| 9 |
} minim_values; |
| 10 |
|
| 11 |
void echarge(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu); |
| 12 |
void echarge1(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu,double **deqq); |
| 13 |
void echarge2(int i,int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,float **hessx,float **hessy,float **hessz); |
| 14 |
|
| 15 |
void echarge(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu) |
| 16 |
{ |
| 17 |
int i, j; |
| 18 |
double xi, yi, zi, xk,yk,zk; |
| 19 |
double xr, yr, zr; |
| 20 |
double e, rik, rik2,cc; |
| 21 |
double rdielc; |
| 22 |
double cutoff; |
| 23 |
|
| 24 |
cutoff = chrgcut*chrgcut; |
| 25 |
*eu = 0.0; |
| 26 |
cc = 4.80298*4.80298*14.39418; |
| 27 |
rdielc = 1.0/units.dielec; |
| 28 |
|
| 29 |
if (minim_values.iprint) |
| 30 |
{ |
| 31 |
fprintf(pcmlogfile,"\nCharge-Charge Interactions : Diele = %7.3f\n",units.dielec); |
| 32 |
fprintf(pcmlogfile," At1 At2 q1 q2 Rik Eqq\n"); |
| 33 |
} |
| 34 |
for (i=1; i < natom; i++) |
| 35 |
{ |
| 36 |
if (charge[i] != 0.0 ) |
| 37 |
{ |
| 38 |
xi = x[i]; |
| 39 |
yi = y[i]; |
| 40 |
zi = z[i]; |
| 41 |
|
| 42 |
for(j=i+1; j <= natom; j++) |
| 43 |
{ |
| 44 |
if (use[i] || use[j]) |
| 45 |
{ |
| 46 |
if ( skip[i][j] != i) |
| 47 |
{ |
| 48 |
if (charge[j] != 0.0) |
| 49 |
{ |
| 50 |
|
| 51 |
xk = x[j]; |
| 52 |
yk = y[j]; |
| 53 |
zk = z[j]; |
| 54 |
|
| 55 |
xr = xi - xk; |
| 56 |
yr = yi - yk; |
| 57 |
zr = zi - zk; |
| 58 |
rik2 = xr*xr + yr*yr + zr*zr; |
| 59 |
if (rik2 < cutoff) |
| 60 |
{ |
| 61 |
rik = sqrt(rik2); |
| 62 |
e = cc*rdielc*charge[i]*charge[j]/rik; |
| 63 |
if (skip[i][j] == -i) |
| 64 |
e /= units.chgscale; |
| 65 |
*eu += e; |
| 66 |
|
| 67 |
if (minim_values.iprint) |
| 68 |
fprintf(pcmlogfile,"QQ: (%-3d)- (%-3d) %-8.3f %-8.3f %-8.3f = %8.4f\n", i, j,charge[i], charge[j],rik,e); |
| 69 |
} |
| 70 |
} |
| 71 |
} |
| 72 |
} |
| 73 |
} |
| 74 |
} |
| 75 |
} |
| 76 |
} |
| 77 |
// ==================================================== |
| 78 |
void echarge1(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu,double **deqq) |
| 79 |
{ |
| 80 |
int i,j; |
| 81 |
double xi, yi, zi; |
| 82 |
double xr, yr, zr; |
| 83 |
double xk,yk,zk; |
| 84 |
double e, rik, cc, rik2; |
| 85 |
double de, dedx, dedy, dedz; |
| 86 |
double rdielc; |
| 87 |
double cutoff; |
| 88 |
|
| 89 |
cutoff = chrgcut*chrgcut; |
| 90 |
cc = 4.80298*4.80298*14.39418; |
| 91 |
rdielc = 1.0/units.dielec; |
| 92 |
|
| 93 |
for (i=1; i < natom; i++) |
| 94 |
{ |
| 95 |
if (charge[i] != 0.0 ) |
| 96 |
{ |
| 97 |
xi = x[i]; |
| 98 |
yi = y[i]; |
| 99 |
zi = z[i]; |
| 100 |
|
| 101 |
for(j=i+1; j <= natom; j++) |
| 102 |
{ |
| 103 |
if (use[i] || use[j]) |
| 104 |
{ |
| 105 |
if ( skip[i][j] != i) |
| 106 |
{ |
| 107 |
if (charge[j] != 0.0) |
| 108 |
{ |
| 109 |
|
| 110 |
xk = x[j]; |
| 111 |
yk = y[j]; |
| 112 |
zk = z[j]; |
| 113 |
xr = xi - xk; |
| 114 |
yr = yi - yk; |
| 115 |
zr = zi - zk; |
| 116 |
rik2 = xr*xr +yr*yr + zr*zr; |
| 117 |
if (rik2 < cutoff) |
| 118 |
{ |
| 119 |
rik = sqrt( rik2); |
| 120 |
e = cc*rdielc*charge[i]*charge[j]/rik; |
| 121 |
if (skip[i][j] == -i) |
| 122 |
e /= units.chgscale; |
| 123 |
|
| 124 |
de = -cc*rdielc*charge[i]*charge[j]/rik2; |
| 125 |
|
| 126 |
de = de/rik; |
| 127 |
dedx = de*xr; |
| 128 |
dedy = de*yr; |
| 129 |
dedz = de*zr; |
| 130 |
|
| 131 |
deqq[i][0] += dedx; |
| 132 |
deqq[i][1] += dedy; |
| 133 |
deqq[i][2] += dedz; |
| 134 |
|
| 135 |
deqq[j][0] -= dedx; |
| 136 |
deqq[j][1] -= dedy; |
| 137 |
deqq[j][2] -= dedz; |
| 138 |
*eu += e; |
| 139 |
} |
| 140 |
} |
| 141 |
} |
| 142 |
} |
| 143 |
} |
| 144 |
} |
| 145 |
} |
| 146 |
} |
| 147 |
// ========================================== |
| 148 |
void echarge2(int i,int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,float **hessx,float **hessy,float **hessz) |
| 149 |
{ |
| 150 |
int j,k; |
| 151 |
double fik,de,d2e,d2edx,d2edy,d2edz; |
| 152 |
double xi,yi,zi,xr,yr,zr,term[3][3]; |
| 153 |
double xk,yk,zk; |
| 154 |
double r,r2; |
| 155 |
double cc, rdielc; |
| 156 |
double cutoff; |
| 157 |
|
| 158 |
if (charge[i] == 0.0) |
| 159 |
return; |
| 160 |
|
| 161 |
cutoff = chrgcut*chrgcut; |
| 162 |
|
| 163 |
cc = 4.80298*4.80298*14.39418; |
| 164 |
rdielc = 1.0/units.dielec; |
| 165 |
|
| 166 |
skip[i][i] = i; |
| 167 |
|
| 168 |
xi = x[i]; |
| 169 |
yi = y[i]; |
| 170 |
zi = z[i]; |
| 171 |
|
| 172 |
|
| 173 |
for (k=1; k <= natom; k++) |
| 174 |
{ |
| 175 |
if (i != k && charge[k] != 0.0 && (skip[i][k] != i) ) |
| 176 |
{ |
| 177 |
xk = x[k]; |
| 178 |
yk = y[k]; |
| 179 |
zk = z[k]; |
| 180 |
|
| 181 |
xr = xi - xk; |
| 182 |
yr = yi - yk; |
| 183 |
zr = zi - zk; |
| 184 |
r2 = xr*xr + yr*yr + zr*zr; |
| 185 |
if (r2 < cutoff) |
| 186 |
{ |
| 187 |
r = sqrt(r2); |
| 188 |
fik = cc*rdielc*charge[i]*charge[k]; |
| 189 |
if (skip[i][k] == -i) |
| 190 |
fik /= units.chgscale; |
| 191 |
de = -fik/r2; |
| 192 |
d2e = -2.0*de/r; |
| 193 |
de = de / r; |
| 194 |
d2e = (d2e-de) / r2; |
| 195 |
d2edx = d2e * xr; |
| 196 |
d2edy = d2e * yr; |
| 197 |
d2edz = d2e * zr; |
| 198 |
term[0][0]= d2edx*xr + de; |
| 199 |
term[1][0] = d2edx*yr; |
| 200 |
term[2][0] = d2edx*zr; |
| 201 |
term[0][1] = term[1][0]; |
| 202 |
term[1][1] = d2edy*yr + de; |
| 203 |
term[2][1] = d2edy*zr; |
| 204 |
term[0][2] = term[2][0]; |
| 205 |
term[1][2] = term[2][1]; |
| 206 |
term[2][2] = d2edz*zr + de; |
| 207 |
|
| 208 |
for (j=0; j < 3; j++) |
| 209 |
{ |
| 210 |
hessx[i][j] += term[j][0]; |
| 211 |
hessy[i][j] += term[j][1]; |
| 212 |
hessz[i][j] += term[j][2]; |
| 213 |
hessx[k][j] -= term[j][0]; |
| 214 |
hessy[k][j] -= term[j][1]; |
| 215 |
hessz[k][j] -= term[j][2]; |
| 216 |
} |
| 217 |
} |
| 218 |
} |
| 219 |
} |
| 220 |
} |
