FastPBC.C

00001 /***************************************************************************
00002 *cr
00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the
00004 *cr University of Illinois
00005 *cr All Rights Reserved
00006 *cr
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010 * RCS INFORMATION:
00011 *
00012 * $RCSfile: FastPBC.C,v $
00013 * $Author: johns $ $Locker: $ $State: Exp $
00014 * $Revision: 1.6 $ $Date: 2022年03月17日 20:21:42 $
00015 *
00016 ***************************************************************************
00017 * DESCRIPTION:
00018 * Fast PBC wrapping code
00019 ***************************************************************************/
00020 #include <stdio.h>
00021 #include <math.h>
00022 #include "Measure.h"
00023 #include "FastPBC.h"
00024 
00025 #ifdef VMDOPENACC
00026 #include <accelmath.h>
00027 #endif
00028 
00033 
00034 // If there is no CUDA defined, fall back to the CPU version.
00035 #ifndef VMDCUDA
00036 void fpbc_exec_unwrap(Molecule* mol, int first, int last, int sellen, int* indexlist) {
00037 fpbc_exec_unwrap_cpu(mol, first, last, sellen, indexlist);
00038 }
00039 
00040 void fpbc_exec_wrapatomic(Molecule* mol, int first, int last, int sellen, int* indexlist, 
00041 float* weights, AtomSel* csel, float* center) {
00042 fpbc_exec_wrapatomic_cpu(mol, first, last, sellen, indexlist, weights, csel, center);
00043 }
00044 
00045 void fpbc_exec_wrapcompound(Molecule* mol, int first, int last, int fnum, int *compoundmap, int sellen, int* indexlist,
00046 float* weights, AtomSel* csel, float* center, float* massarr) {
00047 fpbc_exec_wrapcompound_cpu(mol, first, last, fnum, compoundmap, sellen, indexlist, weights, csel, center, massarr);
00048 }
00049 
00050 void fpbc_exec_join(Molecule* mol, int first, int last, int fnum, int *compoundmap, int sellen, int* indexlist) {
00051 fpbc_exec_join_cpu(mol, first, last, fnum, compoundmap, sellen, indexlist);
00052 }
00053 
00054 void fpbc_exec_recenter(Molecule* mol, int first, int last, int csellen, int* cindexlist, int fnum, int *compoundmap, int sellen, int* indexlist, float* weights, AtomSel* csel, float* massarr) {
00055 fpbc_exec_recenter_cpu(mol, first, last, csellen, cindexlist, fnum, compoundmap, sellen, indexlist, weights, csel, massarr);
00056 }
00057 #endif
00058 
00059 void fpbc_exec_wrapatomic_cpu(Molecule* mol, int first, int last, int sellen, int* indexlist, 
00060 float* weights, AtomSel* csel, float* center) {
00061 int f, i, j, k;
00062 float boxsize[3];
00063 float invboxsize[3];
00064 for (f=first; f<=last; f++) {
00065 Timestep *ts = mol->get_frame(f);
00066 boxsize[0] = ts->a_length;
00067 boxsize[1] = ts->b_length;
00068 boxsize[2] = ts->c_length;
00069 for (j=0; j < 3; j++) {
00070 invboxsize[j] = 1.0f/boxsize[j];
00071 }
00072 //If the center of mass needs to be found, find it! There is a function for that. :D
00073 if (csel != NULL) {
00074 measure_center(csel, ts->pos, weights, center);
00075 }
00076 for (k=0; k<sellen; k++) {
00077 i = indexlist[k];
00078 for (j=0; j < 3; j++) {
00079 //Compute the shift in terms of the number of box-lengths, scale by it, and reposition.
00080 ts->pos[i*3+j] = ts->pos[i*3+j] - (roundf((ts->pos[i*3+j] - center[j]) * invboxsize[j]) * boxsize[j]);
00081 }
00082 }
00083 }
00084 }
00085 
00086 
00087 void fpbc_exec_wrapcompound_cpu(Molecule* mol, int first, int last, int fnum, int *compoundmap, int sellen, int* indexlist,
00088 float* weights, AtomSel* csel, float* center, float* massarr) {
00089 int f, i, j, k, l;
00090 float boxsize[3];
00091 float invboxsize[3];
00092 for (f=first; f<=last; f++) {
00093 Timestep *ts = mol->get_frame(f);
00094 boxsize[0] = ts->a_length;
00095 boxsize[1] = ts->b_length;
00096 boxsize[2] = ts->c_length;
00097 for (j=0; j < 3; j++) {
00098 invboxsize[j] = 1.0f/boxsize[j];
00099 }
00100 //If the center of mass needs to be found, find it! There is a function for that. :D
00101 if (csel != NULL) {
00102 measure_center(csel, ts->pos, weights, center);
00103 }
00104 for (l = 0; l < fnum; l++) {
00105 int lowbound = compoundmap[l];
00106 int highbound = compoundmap[l+1];
00107 //calculate the mass weighted center of the compound
00108 float cmass = 0; // initialize mass accumulator
00109 float ccenter[3] = {0,0,0}; // initialize position accumulator
00110 // cycle through atoms
00111 for (k = lowbound; k < highbound; k++ ) {
00112 i = indexlist[k];
00113 for (j=0; j < 3; j++) {
00114 ccenter[j] += massarr[i] * (ts->pos[i*3+j]);
00115 }
00116 cmass += massarr[i];
00117 }
00118 cmass = 1.0f / cmass;
00119 // divide pos by mass to get final mass weighted com
00120 for (j=0; j < 3; j++) {
00121 ccenter[j] *= cmass;
00122 }
00123 //move the compound
00124 for (k = lowbound; k < highbound; k++ ) {
00125 i = indexlist[k];
00126 for (j=0; j < 3; j++) {
00127 ts->pos[i*3+j] = ts->pos[i*3+j] - (roundf((ccenter[j] - center[j]) * invboxsize[j]) * boxsize[j]);
00128 }
00129 }
00130 }
00131 }
00132 }
00133 
00134 
00135 void fpbc_exec_join_cpu(Molecule* mol, int first, int last, int fnum, int *compoundmap, int sellen, int* indexlist) {
00136 int f, i, j, k, l;
00137 float boxsize[3];
00138 float invboxsize[3];
00139 float *pos;
00140 Timestep *ts;
00141 #pragma acc data copyin (indexlist[sellen], compoundmap[fnum+1])
00142 {
00143 //Loop over all the frames
00144 for (f=first; f<=last; f++) {
00145 //Grab the current coordinates and boxlength.
00146 ts = mol->get_frame(f);
00147 boxsize[0] = ts->a_length;
00148 boxsize[1] = ts->b_length;
00149 boxsize[2] = ts->c_length;
00150 pos = ts->pos;
00151 #pragma acc data copy(pos[3*mol->nAtoms]) copyin(boxsize[3]) create(invboxsize[3]) async(f%2)
00152 {
00153 #pragma acc kernels private(i, j, k, l) async(f%2)
00154 {
00155 //Minimize divisions by computing an inverse box-length.
00156 for (int q=0; q < 3; q++) {
00157 invboxsize[q] = 1.0f/boxsize[q];
00158 }
00159 #pragma acc loop independent private(i, j, k, l)
00160 for (l = 0; l < fnum; l++) {
00161 float ccenter[3];
00162 int lowbound = compoundmap[l];
00163 int highbound = compoundmap[l+1];
00164 i = indexlist[lowbound];
00165 //Use the first element within the compound as the center.
00166 for (j=0; j < 3; j++) {
00167 ccenter[j] = pos[(i+((highbound-lowbound)/2))*3+j];
00168 }
00169 //move the compound, wrapping it to be within half a box dimension from the center
00170 for (k = lowbound; k < highbound; k++ ) {
00171 i = indexlist[k];
00172 for (j=0; j < 3; j++) {
00173 pos[i*3+j] = pos[i*3+j] - (rintf((pos[i*3+j] - ccenter[j]) * invboxsize[j]) * boxsize[j]);
00174 }
00175 }
00176 }
00177 }
00178 }
00179 }
00180 #pragma acc wait
00181 }
00182 }
00183 
00184 
00185 void fpbc_exec_unwrap_cpu(Molecule* mol, int first, int last, int sellen, int* indexlist) {
00186 Timestep *ts, *prev;
00187 int f, i, j, idx;
00188 float boxsize[3];
00189 float oldboxsize[3];
00190 float invboxsize[3];
00191 float *prevw = new float[3*sellen];//previous wrapped coordinates
00192 float curw; //holder for current wrapped coordinates
00193 float disp;//displacement holder
00194 
00195 //Setup previous wrapped coordinates for first go round.
00196 ts = mol->get_frame(first);
00197 for (i=0; i<sellen; i++) {
00198 for (j=0; j<3; j++) {
00199 prevw[3*i+j] = ts->pos[indexlist[i]*3+j];
00200 }
00201 }
00202 
00203 for (f=first+1; f<=last; f++) {
00204 ts = mol->get_frame(f);//current coordinates (wrapped)
00205 prev = mol->get_frame(f-1);//previous unwrapped coordinates
00206 boxsize[0] = ts->a_length;
00207 boxsize[1] = ts->b_length;
00208 boxsize[2] = ts->c_length;
00209 oldboxsize[0] = prev->a_length;
00210 oldboxsize[1] = prev->b_length;
00211 oldboxsize[2] = prev->c_length;
00212 for (j=0; j < 3; j++) {
00213 invboxsize[j] = 1.0f/boxsize[j];
00214 }
00215 for (i=0; i<sellen; i++) {
00216 idx = indexlist[i];
00217 for (j=0; j < 3; j++) {
00218 disp = ts->pos[idx*3+j] - prevw[3*i+j];
00219 curw = ts->pos[idx*3+j];
00220 ts->pos[idx*3+j] = prev->pos[idx*3+j] + disp
00221 -(floorf(disp * invboxsize[j] + 0.5f) * boxsize[j])
00222 -(floorf(((prevw[3*i+j]-prev->pos[idx*3+j])/oldboxsize[j])+0.5f)*(boxsize[j]-oldboxsize[j]));
00223 prevw[3*i+j] = curw;
00224 }
00225 }
00226 }
00227 delete [] prevw;
00228 }
00229 
00230 
00231 void fpbc_exec_recenter_cpu(Molecule* mol, int first, int last, int csellen, int* cindexlist, int fnum, int *compoundmap, int sellen, int* indexlist, float* weights, AtomSel* csel, float* massarr) {
00232 float center[3];//This needs to be passed, but since csel is guaranteed to point somewhere,
00233 //it does not need to be set.
00234 fpbc_exec_unwrap_cpu(mol, first, last, csellen, cindexlist);
00235 if (fnum)
00236 fpbc_exec_wrapcompound_cpu(mol, first, last, fnum, compoundmap, sellen, indexlist, weights, csel, &center[0], massarr);
00237 else
00238 fpbc_exec_wrapatomic_cpu(mol, first, last, sellen, indexlist, weights, csel, &center[0]);
00239 }
00240 
00241 
00242 
00243