|
Lines 1-5
Link Here
|
| 1 |
/* |
1 |
/* |
| 2 |
* $Id: pme.c,v 1.61.2.4 2005/10/06 09:58:08 lindahl Exp $ |
2 |
* $Id: pme.c,v 1.61.2.5 2005/10/13 18:38:04 lindahl Exp $ |
| 3 |
* |
3 |
* |
| 4 |
* This source code is part of |
4 |
* This source code is part of |
| 5 |
* |
5 |
* |
|
Lines 521-528
Link Here
|
| 521 |
range_check(zidx,0,nz); |
521 |
range_check(zidx,0,nz); |
| 522 |
#endif |
522 |
#endif |
| 523 |
i0 = ii0+xidx; /* Pointer arithmetic */ |
523 |
i0 = ii0+xidx; /* Pointer arithmetic */ |
| 524 |
norder = n*4; |
524 |
norder = n*order; |
| 525 |
norder1 = norder+4; |
525 |
norder1 = norder+order; |
| 526 |
|
526 |
|
| 527 |
i = ii0[xidx]; |
527 |
i = ii0[xidx]; |
| 528 |
j = jj0[yidx]; |
528 |
j = jj0[yidx]; |
|
Lines 922-1007
Link Here
|
| 922 |
int nz,rvec fractx[],ivec idx[],real charge[],int nr) |
922 |
int nz,rvec fractx[],ivec idx[],real charge[],int nr) |
| 923 |
{ |
923 |
{ |
| 924 |
/* construct splines for local atoms */ |
924 |
/* construct splines for local atoms */ |
| 925 |
int i,k,l; |
925 |
int i,j,k,l; |
| 926 |
real drXX,drYY,drZZ; |
926 |
real drXX,drYY,drZZ; |
| 927 |
real div,rcons; |
927 |
real dr,div,rcons; |
| 928 |
real *dataXX,*dataYY,*dataZZ; |
928 |
real *dataXX,*dataYY,*dataZZ; |
| 929 |
real *ddataXX,*ddataYY,*ddataZZ; |
929 |
real *ddataXX,*ddataYY,*ddataZZ; |
| 930 |
real tmpX,tmpY,tmpZ,lastX,lastY,lastZ; |
930 |
real tmpX,tmpY,tmpZ,lastX,lastY,lastZ; |
| 931 |
|
931 |
real *data,*ddata,*xptr; |
| 932 |
for(i=0; (i<nr); i++) { |
|
|
| 933 |
|
932 |
|
| 934 |
if (charge[i] != 0.0) { |
933 |
if( order == 4) |
|
|
934 |
{ |
| 935 |
for(i=0; (i<nr); i++) |
| 936 |
{ |
| 935 |
|
937 |
|
| 936 |
drXX = fractx[i][XX]; |
938 |
if (charge[i] != 0.0) |
| 937 |
drYY = fractx[i][YY]; |
939 |
{ |
| 938 |
drZZ = fractx[i][ZZ]; |
|
|
| 939 |
|
| 940 |
/* dr is relative offset from lower cell limit */ |
| 941 |
dataXX=theta[XX]+i*4; |
| 942 |
dataYY=theta[YY]+i*4; |
| 943 |
dataZZ=theta[ZZ]+i*4; |
| 944 |
|
| 945 |
dataXX[3]=0; |
| 946 |
dataYY[3]=0; |
| 947 |
dataZZ[3]=0; |
| 948 |
dataXX[1]=drXX; |
| 949 |
dataYY[1]=drYY; |
| 950 |
dataZZ[1]=drZZ; |
| 951 |
dataXX[0]=1.0-drXX; |
| 952 |
dataYY[0]=1.0-drYY; |
| 953 |
dataZZ[0]=1.0-drZZ; |
| 954 |
|
| 955 |
dataXX[2]=0.5*drXX*dataXX[1]; |
| 956 |
dataYY[2]=0.5*drYY*dataYY[1]; |
| 957 |
dataZZ[2]=0.5*drZZ*dataZZ[1]; |
| 958 |
|
| 959 |
dataXX[1]=0.5*((drXX+1.0)*dataXX[0]+(2.0-drXX)*dataXX[1]); |
| 960 |
dataYY[1]=0.5*((drYY+1.0)*dataYY[0]+(2.0-drYY)*dataYY[1]); |
| 961 |
dataZZ[1]=0.5*((drZZ+1.0)*dataZZ[0]+(2.0-drZZ)*dataZZ[1]); |
| 962 |
|
| 963 |
dataXX[0]=0.5*(1.0-drXX)*dataXX[0]; |
| 964 |
dataYY[0]=0.5*(1.0-drYY)*dataYY[0]; |
| 965 |
dataZZ[0]=0.5*(1.0-drZZ)*dataZZ[0]; |
| 966 |
|
| 967 |
/* differentiate */ |
| 968 |
ddataXX = dtheta[XX]+i*4; |
| 969 |
ddataYY = dtheta[YY]+i*4; |
| 970 |
ddataZZ = dtheta[ZZ]+i*4; |
| 971 |
|
| 972 |
ddataXX[0] = -dataXX[0]; |
| 973 |
ddataYY[0] = -dataYY[0]; |
| 974 |
ddataZZ[0] = -dataZZ[0]; |
| 975 |
ddataXX[1] = dataXX[0]-dataXX[1]; |
| 976 |
ddataYY[1] = dataYY[0]-dataYY[1]; |
| 977 |
ddataZZ[1] = dataZZ[0]-dataZZ[1]; |
| 978 |
ddataXX[2] = dataXX[1]-dataXX[2]; |
| 979 |
ddataYY[2] = dataYY[1]-dataYY[2]; |
| 980 |
ddataZZ[2] = dataZZ[1]-dataZZ[2]; |
| 981 |
ddataXX[3] = dataXX[2]-dataXX[3]; |
| 982 |
ddataYY[3] = dataYY[2]-dataYY[3]; |
| 983 |
ddataZZ[3] = dataZZ[2]-dataZZ[3]; |
| 984 |
|
| 985 |
div=1.0/3.0; |
| 986 |
dataXX[3]=div*drXX*dataXX[2]; |
| 987 |
dataYY[3]=div*drYY*dataYY[2]; |
| 988 |
dataZZ[3]=div*drZZ*dataZZ[2]; |
| 989 |
|
| 990 |
dataXX[2]=div*((drXX+1.0)*dataXX[1]+(3.0-drXX)*dataXX[2]); |
| 991 |
dataYY[2]=div*((drYY+1.0)*dataYY[1]+(3.0-drYY)*dataYY[2]); |
| 992 |
dataZZ[2]=div*((drZZ+1.0)*dataZZ[1]+(3.0-drZZ)*dataZZ[2]); |
| 993 |
|
| 994 |
dataXX[1]=div*((drXX+2.0)*dataXX[0]+(2.0-drXX)*dataXX[1]); |
| 995 |
dataYY[1]=div*((drYY+2.0)*dataYY[0]+(2.0-drYY)*dataYY[1]); |
| 996 |
dataZZ[1]=div*((drZZ+2.0)*dataZZ[0]+(2.0-drZZ)*dataZZ[1]); |
| 997 |
|
| 998 |
dataXX[0]=div*(1.0-drXX)*dataXX[0]; |
| 999 |
dataYY[0]=div*(1.0-drYY)*dataYY[0]; |
| 1000 |
dataZZ[0]=div*(1.0-drZZ)*dataZZ[0]; |
| 1001 |
|
940 |
|
|
|
941 |
drXX = fractx[i][XX]; |
| 942 |
drYY = fractx[i][YY]; |
| 943 |
drZZ = fractx[i][ZZ]; |
| 944 |
|
| 945 |
/* dr is relative offset from lower cell limit */ |
| 946 |
dataXX=theta[XX]+i*4; |
| 947 |
dataYY=theta[YY]+i*4; |
| 948 |
dataZZ=theta[ZZ]+i*4; |
| 949 |
|
| 950 |
dataXX[3]=0; |
| 951 |
dataYY[3]=0; |
| 952 |
dataZZ[3]=0; |
| 953 |
dataXX[1]=drXX; |
| 954 |
dataYY[1]=drYY; |
| 955 |
dataZZ[1]=drZZ; |
| 956 |
dataXX[0]=1.0-drXX; |
| 957 |
dataYY[0]=1.0-drYY; |
| 958 |
dataZZ[0]=1.0-drZZ; |
| 959 |
|
| 960 |
dataXX[2]=0.5*drXX*dataXX[1]; |
| 961 |
dataYY[2]=0.5*drYY*dataYY[1]; |
| 962 |
dataZZ[2]=0.5*drZZ*dataZZ[1]; |
| 963 |
|
| 964 |
dataXX[1]=0.5*((drXX+1.0)*dataXX[0]+(2.0-drXX)*dataXX[1]); |
| 965 |
dataYY[1]=0.5*((drYY+1.0)*dataYY[0]+(2.0-drYY)*dataYY[1]); |
| 966 |
dataZZ[1]=0.5*((drZZ+1.0)*dataZZ[0]+(2.0-drZZ)*dataZZ[1]); |
| 967 |
|
| 968 |
dataXX[0]=0.5*(1.0-drXX)*dataXX[0]; |
| 969 |
dataYY[0]=0.5*(1.0-drYY)*dataYY[0]; |
| 970 |
dataZZ[0]=0.5*(1.0-drZZ)*dataZZ[0]; |
| 971 |
|
| 972 |
/* differentiate */ |
| 973 |
ddataXX = dtheta[XX]+i*4; |
| 974 |
ddataYY = dtheta[YY]+i*4; |
| 975 |
ddataZZ = dtheta[ZZ]+i*4; |
| 976 |
|
| 977 |
ddataXX[0] = -dataXX[0]; |
| 978 |
ddataYY[0] = -dataYY[0]; |
| 979 |
ddataZZ[0] = -dataZZ[0]; |
| 980 |
ddataXX[1] = dataXX[0]-dataXX[1]; |
| 981 |
ddataYY[1] = dataYY[0]-dataYY[1]; |
| 982 |
ddataZZ[1] = dataZZ[0]-dataZZ[1]; |
| 983 |
ddataXX[2] = dataXX[1]-dataXX[2]; |
| 984 |
ddataYY[2] = dataYY[1]-dataYY[2]; |
| 985 |
ddataZZ[2] = dataZZ[1]-dataZZ[2]; |
| 986 |
ddataXX[3] = dataXX[2]-dataXX[3]; |
| 987 |
ddataYY[3] = dataYY[2]-dataYY[3]; |
| 988 |
ddataZZ[3] = dataZZ[2]-dataZZ[3]; |
| 989 |
|
| 990 |
div=1.0/3.0; |
| 991 |
dataXX[3]=div*drXX*dataXX[2]; |
| 992 |
dataYY[3]=div*drYY*dataYY[2]; |
| 993 |
dataZZ[3]=div*drZZ*dataZZ[2]; |
| 994 |
|
| 995 |
dataXX[2]=div*((drXX+1.0)*dataXX[1]+(3.0-drXX)*dataXX[2]); |
| 996 |
dataYY[2]=div*((drYY+1.0)*dataYY[1]+(3.0-drYY)*dataYY[2]); |
| 997 |
dataZZ[2]=div*((drZZ+1.0)*dataZZ[1]+(3.0-drZZ)*dataZZ[2]); |
| 998 |
|
| 999 |
dataXX[1]=div*((drXX+2.0)*dataXX[0]+(2.0-drXX)*dataXX[1]); |
| 1000 |
dataYY[1]=div*((drYY+2.0)*dataYY[0]+(2.0-drYY)*dataYY[1]); |
| 1001 |
dataZZ[1]=div*((drZZ+2.0)*dataZZ[0]+(2.0-drZZ)*dataZZ[1]); |
| 1002 |
|
| 1003 |
dataXX[0]=div*(1.0-drXX)*dataXX[0]; |
| 1004 |
dataYY[0]=div*(1.0-drYY)*dataYY[0]; |
| 1005 |
dataZZ[0]=div*(1.0-drZZ)*dataZZ[0]; |
| 1006 |
} |
| 1002 |
} |
1007 |
} |
| 1003 |
} |
1008 |
} |
| 1004 |
|
1009 |
else |
|
|
1010 |
{ |
| 1011 |
/* general case, order != 4 */ |
| 1012 |
for(i=0; (i<nr); i++) |
| 1013 |
{ |
| 1014 |
if (charge[i] != 0.0) |
| 1015 |
{ |
| 1016 |
xptr = fractx[i]; |
| 1017 |
for(j=0; (j<DIM); j++) |
| 1018 |
{ |
| 1019 |
dr = xptr[j]; |
| 1020 |
|
| 1021 |
/* dr is relative offset from lower cell limit */ |
| 1022 |
data=&(theta[j][i*order]); |
| 1023 |
data[order-1]=0; |
| 1024 |
data[1]=dr; |
| 1025 |
data[0]=1-dr; |
| 1026 |
|
| 1027 |
for(k=3; (k<order); k++) |
| 1028 |
{ |
| 1029 |
div=1.0/(k-1.0); |
| 1030 |
data[k-1]=div*dr*data[k-2]; |
| 1031 |
for(l=1; (l<(k-1)); l++) |
| 1032 |
{ |
| 1033 |
data[k-l-1]=div*((dr+l)*data[k-l-2]+(k-l-dr)*data[k-l-1]); |
| 1034 |
} |
| 1035 |
data[0]=div*(1-dr)*data[0]; |
| 1036 |
} |
| 1037 |
/* differentiate */ |
| 1038 |
ddata = &(dtheta[j][i*order]); |
| 1039 |
ddata[0] = -data[0]; |
| 1040 |
for(k=1; (k<order); k++) |
| 1041 |
{ |
| 1042 |
ddata[k]=data[k-1]-data[k]; |
| 1043 |
} |
| 1044 |
div=1.0/(order-1); |
| 1045 |
data[order-1]=div*dr*data[order-2]; |
| 1046 |
for(l=1; (l<(order-1)); l++) |
| 1047 |
{ |
| 1048 |
data[order-l-1]=div*((dr+l)*data[order-l-2]+(order-l-dr)*data[order-l-1]); |
| 1049 |
} |
| 1050 |
data[0]=div*(1-dr)*data[0]; |
| 1051 |
} |
| 1052 |
} |
| 1053 |
} |
| 1054 |
} |
| 1005 |
} |
1055 |
} |
| 1006 |
|
1056 |
|
| 1007 |
|
1057 |
|
