00001
00007 #include <string.h>
00008 #include "OptPmeRealSpace.h"
00009 #include <assert.h>
00010
00011 OptPmeRealSpace::OptPmeRealSpace(PmeGrid grid,int natoms)
00012 : myGrid(grid), N(natoms) {
00013 int order = myGrid.order;
00014 M = new double[3*N*order];
00015 dM = new double[3*N*order];
00016 }
00017
00018 OptPmeRealSpace::~OptPmeRealSpace() {
00019 delete [] M;
00020 delete [] dM;
00021 }
00022
00023
00024
00025 #define order 4
00026
00027 void OptPmeRealSpace::fill_charges(double **q_arr, PmeParticle p[], int zstart, int zlen) {
00028
00029 int i, j, k, l;
00030 int stride;
00031 int K1, K2, K3, dim2, dim3;
00032 int32 K3_1;
00033 double *Mi, *dMi;
00034 Mi = M; dMi = dM;
00035 K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2; dim3=myGrid.dim3;
00036 K3_1 = K3 - 1;
00037
00038 stride = 3*order;
00039
00040 for (i=0; i<N; i++) {
00041 double q;
00042 int32 u1, u2, u2i, u3i;
00043 double fr1 = p[i].x;
00044 double fr2 = p[i].y;
00045 double fr3 = p[i].z;
00046 q = p[i].cg;
00047 u1 = (int)fr1;
00048 u2i = (int)fr2;
00049 u3i = (int)fr3;
00050 fr1 -= u1;
00051 fr2 -= u2i;
00052 fr3 -= u3i;
00053
00054
00055 Mi[0] = ( ( (-1./6.) * fr1 + 0.5 ) * fr1 - 0.5 ) * fr1 + (1./6.);
00056 Mi[1] = ( ( 0.5 * fr1 - 1.0 ) * fr1 ) * fr1 + (2./3.);
00057 Mi[2] = ( ( -0.5 * fr1 + 0.5 ) * fr1 + 0.5 ) * fr1 + (1./6.);
00058 Mi[3] = (1./6.) * fr1 * fr1 * fr1;
00059 dMi[0] = ( -0.5 * fr1 + 1.0 )* fr1 - 0.5;
00060 dMi[1] = ( 1.5 * fr1 - 2.0 ) * fr1;
00061 dMi[2] = ( -1.5 * fr1 + 1.0 ) * fr1 + 0.5;
00062 dMi[3] = 0.5 * fr1 * fr1;
00063 Mi[4] = ( ( (-1./6.) * fr2 + 0.5 ) * fr2 - 0.5 ) * fr2 + (1./6.);
00064 Mi[5] = ( ( 0.5 * fr2 - 1.0 ) * fr2 ) * fr2 + (2./3.);
00065 Mi[6] = ( ( -0.5 * fr2 + 0.5 ) * fr2 + 0.5 ) * fr2 + (1./6.);
00066 Mi[7] = (1./6.) * fr2 * fr2 * fr2;
00067 dMi[4] = ( -0.5 * fr2 + 1.0 )* fr2 - 0.5;
00068 dMi[5] = ( 1.5 * fr2 - 2.0 ) * fr2;
00069 dMi[6] = ( -1.5 * fr2 + 1.0 ) * fr2 + 0.5;
00070 dMi[7] = 0.5 * fr2 * fr2;
00071 Mi[8] = ( ( (-1./6.) * fr3 + 0.5 ) * fr3 - 0.5 ) * fr3 + (1./6.);
00072 Mi[9] = ( ( 0.5 * fr3 - 1.0 ) * fr3 ) * fr3 + (2./3.);
00073 Mi[10] = ( ( -0.5 * fr3 + 0.5 ) * fr3 + 0.5 ) * fr3 + (1./6.);
00074 Mi[11] = (1./6.) * fr3 * fr3 * fr3;
00075 dMi[8] = ( -0.5 * fr3 + 1.0 )* fr3 - 0.5;
00076 dMi[9] = ( 1.5 * fr3 - 2.0 ) * fr3;
00077 dMi[10] = ( -1.5 * fr3 + 1.0 ) * fr3 + 0.5;
00078 dMi[11] = 0.5 * fr3 * fr3;
00079
00080 u1 -= order;
00081 u2i -= order;
00082 u3i -= order;
00083 u3i += 1;
00084 for (j=0; j<order; j++) {
00085 double m1;
00086 int32 ind1;
00087 m1 = Mi[j]*q;
00088 u1++;
00089
00090 ind1 = (u1 + ((unsigned)u1 >>31)*K1)*dim2;
00091 u2 = u2i;
00092 for (k=0; k<order; k++) {
00093 double m1m2;
00094 int32 ind2;
00095 m1m2 = m1*Mi[order+k];
00096 u2++;
00097
00098
00099 ind2 = ind1 + (u2 + ((unsigned)u2 >>31)*K2);
00100
00101 double *qline = q_arr[ind2];
00102
00103
00104
00105
00106
00107 for (l=0; l<order; l++) {
00108 double m3;
00109 int32 ind;
00110 m3 = Mi[2*order + l];
00111 int32 u3 = u3i + l;
00112 ind = u3 - zstart;
00113
00114 ind = ind - ((unsigned)(K3_1 - ind) >>31)*K3;
00115
00116
00117
00118
00119
00120 qline[ind] += m1m2*m3;
00121 }
00122 }
00123 }
00124 Mi += stride;
00125 dMi += stride;
00126 }
00127 }
00128
00129 void OptPmeRealSpace::compute_forces(const double * const *q_arr,
00130 const PmeParticle p[], Vector f[],
00131 int zstart, int zlen) {
00132
00133 int i, j, k, l, stride;
00134 double f1, f2, f3;
00135 double *Mi, *dMi;
00136 int K1, K2, K3, dim2;
00137 int32 K3_1;
00138
00139 K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2;
00140 K3_1 = K3 - 1;
00141
00142 stride=3*order;
00143 Mi = M; dMi = dM;
00144
00145 for (i=0; i<N; i++) {
00146 double q;
00147 int32 u1, u2, u2i, u3i;
00148 q = p[i].cg;
00149 f1=f2=f3=0.0;
00150 u1 = (int)(p[i].x);
00151 u2i = (int)(p[i].y);
00152 u3i = (int)(p[i].z);
00153 u1 -= order;
00154 u2i -= order;
00155 u3i -= order;
00156 u3i += 1;
00157 for (j=0; j<order; j++) {
00158 double m1, d1;
00159 int32 ind1;
00160 m1=Mi[j]*q;
00161 d1=K1*dMi[j]*q;
00162 u1++;
00163
00164 ind1 = (u1 + ((unsigned)u1 >>31)*K1)*dim2;
00165 u2 = u2i;
00166 for (k=0; k<order; k++) {
00167 double m2, d2, m1m2, m1d2, d1m2;
00168 int32 ind2;
00169 m2=Mi[order+k];
00170 d2=K2*dMi[order+k];
00171 m1m2=m1*m2;
00172 m1d2=m1*d2;
00173 d1m2=d1*m2;
00174 u2++;
00175
00176 ind2 = ind1 + (u2 + ((unsigned)u2 >>31)*K2);
00177 const double *qline = q_arr[ind2];
00178 for (l=0; l<order; l++) {
00179 double term, m3, d3;
00180 int32 ind;
00181 m3=Mi[2*order+l];
00182 d3=K3*dMi[2*order+l];
00183 int32 u3 = u3i + l;
00184 ind = u3 - zstart;
00185
00186 ind = ind - ((unsigned)(K3_1 -ind) >>31)*K3;
00187
00188
00189
00190
00191
00192 term = qline[ind];
00193 f1 -= d1m2 * m3 * term;
00194 f2 -= m1d2 * m3 * term;
00195 f3 -= m1m2 * d3 * term;
00196 }
00197 }
00198 }
00199 Mi += stride;
00200 dMi += stride;
00201 f[i].x = f1;
00202 f[i].y = f2;
00203 f[i].z = f3;
00204 }
00205 }