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: Timestep.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.74 $ $Date: 2020年10月21日 05:43:54 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * 00019 * The Timestep class, which stores coordinates, energies, etc. for a 00020 * single timestep. 00021 * 00022 * Note: As more data is stored for each step, it should go in here. For 00023 * example, H-Bonds could be calculated each step. 00024 ***************************************************************************/ 00025 00026 #include <math.h> 00027 #include "Timestep.h" 00028 #include "Inform.h" 00029 #include "utilities.h" 00030 00031 // allocate memory and return a pointer that is aligned on a given 00032 // byte boundary, to be used for page- or sector-aligned I/O buffers 00033 // We use this since posix_memalign() is not widely available... 00034 #if defined(_WIN64) /* sizeof(size_t) == sizeof(void*) */ 00035 #define myintptrtype size_t 00036 #elif 1 /* sizeof(unsigned long) == sizeof(void*) */ 00037 #define myintptrtype unsigned long 00038 #else 00039 #define myintptrtype uintptr_t /* C99 */ 00040 #endif 00041 00042 00043 static void *alloc_aligned_ptr(size_t sz, size_t blocksz, void **unalignedptr) { 00044 // pad the allocation to an even multiple of the block size 00045 size_t padsz = (sz + (blocksz - 1)) & (~(blocksz - 1)); 00046 void * ptr = calloc(1, padsz + blocksz + blocksz); 00047 *unalignedptr = ptr; 00048 return (void *) ((((myintptrtype) ptr) + (blocksz-1)) & (~(blocksz-1))); 00049 } 00050 00051 00053 Timestep::Timestep(int n, int pagealignsz) { 00054 for(int i=0; i < TSENERGIES; energy[i++] = 0.0); 00055 num = n; 00056 page_align_sz = pagealignsz; 00057 00058 if (page_align_sz > 1) { 00059 pos = (float *) alloc_aligned_ptr(3L*num*sizeof(float), 00060 page_align_sz, (void**) &pos_ptr); 00061 } else { 00062 pos_ptr = (float *) calloc(1, 3L*num*sizeof(float)); 00063 pos=pos_ptr; 00064 } 00065 vel = NULL; 00066 force = NULL; 00067 user = NULL; 00068 user2 = NULL; 00069 user3 = NULL; 00070 user4 = NULL; 00071 qm_timestep = NULL; 00072 a_length = b_length = c_length = 0; 00073 alpha = beta = gamma = 90; 00074 timesteps=0; 00075 physical_time=0; 00076 } 00077 00078 00080 Timestep::Timestep(const Timestep& ts) { 00081 num = ts.num; 00082 page_align_sz = ts.page_align_sz; 00083 00084 // If we support block-based direct I/O, we must use memory buffers 00085 // that are padded to a full block size, and we have to track the 00086 // raw pointer for use at deallocation time. 00087 if (page_align_sz > 1) { 00088 pos = (float *) alloc_aligned_ptr(3L*num*sizeof(float), 00089 page_align_sz, (void**) &pos_ptr); 00090 } else { 00091 pos_ptr = (float *) calloc(1, 3L*num*sizeof(float)); 00092 pos=pos_ptr; 00093 } 00094 memcpy(pos, ts.pos, 3L*num*sizeof(float)); // copy only payload data 00095 00096 if (ts.force) { 00097 force = new float[3L * num]; 00098 memcpy(force, ts.force, 3L * num * sizeof(float)); 00099 } else { 00100 force = NULL; 00101 } 00102 00103 if (ts.vel) { 00104 vel = new float[3L * num]; 00105 memcpy(vel, ts.vel, 3L * num * sizeof(float)); 00106 } else { 00107 vel = NULL; 00108 } 00109 00110 if (ts.user) { 00111 user = new float[num]; 00112 memcpy(user, ts.user, num*sizeof(float)); 00113 } else { 00114 user = NULL; 00115 } 00116 00117 if (ts.user2) { 00118 user2 = new float[num]; 00119 memcpy(user2, ts.user2, num*sizeof(float)); 00120 } else { 00121 user2 = NULL; 00122 } 00123 00124 if (ts.user3) { 00125 user3 = new float[num]; 00126 memcpy(user3, ts.user3, num*sizeof(float)); 00127 } else { 00128 user3 = NULL; 00129 } 00130 00131 if (ts.user4) { 00132 user4 = new float[num]; 00133 memcpy(user4, ts.user4, num*sizeof(float)); 00134 } else { 00135 user4 = NULL; 00136 } 00137 00138 if (ts.qm_timestep) { 00139 qm_timestep = new QMTimestep(*(ts.qm_timestep)); 00140 } else { 00141 qm_timestep = NULL; 00142 } 00143 00144 memcpy(energy, ts.energy, sizeof(ts.energy)); 00145 a_length = ts.a_length; 00146 b_length = ts.b_length; 00147 c_length = ts.c_length; 00148 alpha = ts.alpha; 00149 beta = ts.beta; 00150 gamma = ts.gamma; 00151 timesteps=ts.timesteps; 00152 physical_time=ts.physical_time; 00153 } 00154 00155 00157 Timestep::~Timestep() { 00158 delete [] force; 00159 delete [] vel; 00160 if (pos_ptr) 00161 free(pos_ptr); 00162 delete [] user; 00163 delete [] user2; 00164 delete [] user3; 00165 delete [] user4; 00166 if (qm_timestep) 00167 delete qm_timestep; 00168 } 00169 00170 00171 // reset coords and related items to 0 00172 void Timestep::zero_values() { 00173 if (num <= 0) 00174 return; 00175 00176 memset(pos, 0, 3L*num*sizeof(float)); 00177 00178 for(int i=0; i < TSENERGIES; energy[i++] = 0.0); 00179 timesteps=0; 00180 } 00181 00182 00183 void Timestep::get_transform_vectors(float A[3], float B[3], float C[3]) const 00184 { 00185 // notes: a, b, c are side lengths of the unit cell 00186 // alpha = angle between b and c 00187 // beta = angle between a and c 00188 // gamma = angle between a and b 00189 00190 // convert from degrees to radians 00191 double cosBC = cos(DEGTORAD(alpha)); 00192 double cosAC = cos(DEGTORAD(beta)); 00193 double cosAB, sinAB; 00194 sincos(DEGTORAD(gamma), &sinAB, &cosAB); 00195 00196 // A will lie along the positive x axis. 00197 // B will lie in the x-y plane 00198 // The origin will be (0,0,0). 00199 float Ax = (float) (a_length); 00200 float Bx = (float) (b_length * cosAB); 00201 float By = (float) (b_length * sinAB); 00202 00203 float Cx=0, Cy=0, Cz=0; 00204 // If sinAB is zero, then we can't determine C uniquely since it's defined 00205 // in terms of the angle between A and B. 00206 if (sinAB > 0) { 00207 Cx = (float) cosAC; 00208 Cy = (float) ((cosBC - cosAC * cosAB) / sinAB); 00209 Cz = sqrtf(1.0f - Cx*Cx - Cy*Cy); 00210 } 00211 Cx *= c_length; 00212 Cy *= c_length; 00213 Cz *= c_length; 00214 vec_zero(A); A[0] = Ax; 00215 vec_zero(B); B[0] = Bx; B[1] = By; 00216 vec_zero(C); C[0] = Cx; C[1] = Cy; C[2] = Cz; 00217 } 00218 00219 00220 void Timestep::get_transforms(Matrix4 &a, Matrix4 &b, Matrix4 &c) const { 00221 float A[3], B[3], C[3]; 00222 get_transform_vectors(A, B, C); 00223 a.translate(A); 00224 b.translate(B); 00225 c.translate(C); 00226 } 00227 00228 00229 void Timestep::get_transform_from_cell(const int *cell, Matrix4 &mat) const { 00230 float A[3], B[3], C[3]; 00231 get_transform_vectors(A, B, C); 00232 float Ax=A[0]; 00233 float Bx=B[0]; 00234 float By=B[1]; 00235 float Cx=C[0]; 00236 float Cy=C[1]; 00237 float Cz=C[2]; 00238 mat.identity(); 00239 mat.mat[12] = cell[0]*Ax + cell[1]*Bx + cell[2]*Cx; 00240 mat.mat[13] = cell[1]*By + cell[2]*Cy; 00241 mat.mat[14] = cell[2]*Cz; 00242 } 00243