23 #include "ComputeMgr.decl.h" 29 #define MIN_DEBUG_LEVEL 3 32 #define USE_GLOBALMASTER_VIRIAL_KERNEL 1 41 DebugM(3,
"Constructing client\n");
47 isRequestedAllocSize = 0;
49 numGroupsRequested = 0;
54 if ( sp->
colvarsOn ) forceSendEnabled = 1;
55 if ( sp->
IMDon ) forceSendEnabled = 1;
67 forcePtrs =
new Force*[numPatches];
68 atomPtrs =
new FullAtom*[numPatches];
69 for (
int i = 0; i < numPatches; ++i ) { forcePtrs[i] = 0; atomPtrs[i] = 0; }
73 mass_soa =
new float*[numPatches];
74 pos_soa_x =
new double*[numPatches];
75 pos_soa_y =
new double*[numPatches];
76 pos_soa_z =
new double*[numPatches];
77 force_soa_x =
new double*[numPatches];
78 force_soa_y =
new double*[numPatches];
79 force_soa_z =
new double*[numPatches];
80 transform_soa_i =
new int*[numPatches];
81 transform_soa_j =
new int*[numPatches];
82 transform_soa_k =
new int*[numPatches];
83 for (
int i = 0; i < numPatches; ++i ) {
88 force_soa_x[i] = NULL;
89 force_soa_y[i] = NULL;
90 force_soa_z[i] = NULL;
91 transform_soa_i[i] = NULL;
92 transform_soa_j[i] = NULL;
93 transform_soa_k[i] = NULL;
103 transform_soa_i = NULL;
104 transform_soa_j = NULL;
105 transform_soa_k = NULL;
107 gridForcesPtrs =
new ForceList **[numPatches];
108 numGridObjects = numActiveGridObjects = 0;
109 for (
int i = 0; i < numPatches; ++i ) {
110 forcePtrs[i] = NULL; atomPtrs[i] = NULL;
111 gridForcesPtrs[i] = NULL;
117 delete[] isRequested;
120 delete[] gridForcesPtrs;
124 if(mass_soa)
delete [] mass_soa;
125 if(pos_soa_x)
delete [] pos_soa_x;
126 if(pos_soa_y)
delete [] pos_soa_y;
127 if(pos_soa_z)
delete [] pos_soa_z;
128 if(force_soa_x)
delete [] force_soa_x;
129 if(force_soa_y)
delete [] force_soa_y;
130 if(force_soa_z)
delete [] force_soa_z;
131 if(transform_soa_i)
delete [] transform_soa_i;
132 if(transform_soa_j)
delete [] transform_soa_j;
133 if(transform_soa_k)
delete [] transform_soa_k;
137 DebugM(4,
"Receiving configuration (" << newaid.
size() <<
138 " atoms, " << newgdef.
size() <<
" atoms/groups and " <<
139 newgridobjid.
size() <<
" grid objects) on client\n" <<
endi);
143 if ( forceSendEnabled ) {
146 for (a=newaid.
begin(),a_e=newaid.
end(); a!=a_e; ++a) {
147 if ( *a > max ) max = *a;
149 for (a=newgdef.
begin(),a_e=newgdef.
end(); a!=a_e; ++a) {
150 if ( *a > max ) max = *a;
152 endRequested = max+1;
153 if ( endRequested > isRequestedAllocSize ) {
154 delete [] isRequested;
155 isRequestedAllocSize = endRequested+10;
156 isRequested =
new char[isRequestedAllocSize];
157 memset(isRequested, 0, isRequestedAllocSize);
159 for (a=aid.
begin(),a_e=aid.
end(); a!=a_e; ++a) {
162 for (a=gdef.
begin(),a_e=gdef.
end(); a!=a_e; ++a) {
163 if ( *a != -1 ) isRequested[*a] = 0;
176 if (newgridobjid.
size()) configureGridObjects(newgridobjid);
178 if ( forceSendEnabled ) {
180 for (a=aid.
begin(),a_e=aid.
end(); a!=a_e; ++a) {
183 for (a=gdef.
begin(),a_e=gdef.
end(); a!=a_e; ++a) {
184 if ( *a == -1 ) ++newgcount;
186 isRequested[*a] |= 2;
190 std::sort(gpair.
begin(),gpair.
end());
191 numGroupsRequested = newgcount;
193 DebugM(3,
"Done configure on client\n");
196 void ComputeGlobal::deleteGridObjects()
198 if (numGridObjects == 0)
return;
200 for (ap = ap.begin(); ap != ap.end(); ap++) {
201 ForceList **gridForces = gridForcesPtrs[ap->p->getPatchID()];
202 if (gridForces != NULL) {
203 for (
size_t ig = 0; ig < numGridObjects; ig++) {
204 if (gridForces[ig] != NULL) {
205 delete gridForces[ig];
206 gridForces[ig] = NULL;
209 delete [] gridForces;
213 numGridObjects = numActiveGridObjects = 0;
216 void ComputeGlobal::configureGridObjects(
IntList &newgridobjid)
223 numActiveGridObjects = 0;
225 gridObjActive.
resize(numGridObjects);
230 for ( ; goid_i != goid_e; goid_i++) {
231 if ((*goid_i < 0) || (*goid_i >= numGridObjects)) {
232 NAMD_bug(
"Requested illegal gridForceGrid index.");
234 DebugM(3,
"Adding grid with index " << *goid_i <<
" to ComputeGlobal\n");
235 gridObjActive[*goid_i] = 1;
236 numActiveGridObjects++;
240 for (
size_t ig = 0; ig < numGridObjects; ig++) {
241 DebugM(3,
"Grid index " << ig <<
" is active or inactive? " 242 << gridObjActive[ig] <<
"\n" <<
endi);
246 for (ap = ap.begin(); ap != ap.end(); ap++) {
247 gridForcesPtrs[ap->p->getPatchID()] =
new ForceList *[numGridObjects];
248 ForceList **gridForces = gridForcesPtrs[ap->p->getPatchID()];
249 for (
size_t ig = 0; ig < numGridObjects; ig++) {
250 if (gridObjActive[ig]) {
253 gridForces[ig] = NULL;
260 void ComputeGlobal::recvConfig(ComputeGlobalConfigMsg *msg) {
261 DebugM(3,
"Receiving configure on client\n");
262 configure(msg->aid,msg->gdef);
269 DebugM(3,
"Receiving results (" << msg->
aid.
size() <<
" forces, " 270 << msg->
newgdef.
size() <<
" new group atoms) on client thread " << CthGetToken(CthSelf())->serialNo <<
" msg->resendCoordinates " << msg->
resendCoordinates <<
" msg->totalforces " << msg->
totalforces<<
"\n");
273 if ( forceSendActive && ! forceSendEnabled )
NAMD_bug(
"ComputeGlobal::recvResults forceSendActive without forceSendEnabled");
283 Force **f = forcePtrs;
288 for (ap = ap.
begin(); ap != ap.
end(); ap++) {
289 (*ap).r = (*ap).forceBox->open();
291 t[(*ap).patchID] = (*ap).p->getAtomList().begin();
296 mass_soa[pId] = (*ap).p->patchDataSOA.mass;
297 force_soa_x[pId] = (*ap).p->patchDataSOA.f_global_x;
298 force_soa_y[pId] = (*ap).p->patchDataSOA.f_global_y;
299 force_soa_z[pId] = (*ap).p->patchDataSOA.f_global_z;
300 transform_soa_i[pId] = (*ap).p->patchDataSOA.transform_i;
301 transform_soa_j[pId] = (*ap).p->patchDataSOA.transform_j;
302 transform_soa_k[pId] = (*ap).p->patchDataSOA.transform_k;
311 for ( ; a != a_e; ++a, ++f2 ) {
325 for ( ; a != a_e; ++a, ++f2 ) {
326 DebugM(1,
"processing atom "<<(*a)<<
", F="<<(*f2)<<
"...\n");
328 localID = atomMap->
localID(*a);
330 lidx = localID.
index;
331 if ( lpid ==
notUsed || ! f[lpid] )
continue;
336 force_soa_x[lpid][lidx] += f_atom.
x;
337 force_soa_y[lpid][lidx] += f_atom.
y;
338 force_soa_z[lpid][lidx] += f_atom.
z;
339 #ifndef USE_GLOBALMASTER_VIRIAL_KERNEL 340 x_orig.
x = pos_soa_x[lpid][lidx];
341 x_orig.
y = pos_soa_y[lpid][lidx];
342 x_orig.
z = pos_soa_z[lpid][lidx];
343 trans.
i = transform_soa_i[lpid][lidx];
344 trans.
j = transform_soa_j[lpid][lidx];
345 trans.
k = transform_soa_k[lpid][lidx];
348 extVirial +=
outer(f_atom,x_atom);
352 for ( ; a != a_e; ++a, ++f2 ) {
353 DebugM(1,
"processing atom "<<(*a)<<
", F="<<(*f2)<<
"...\n");
357 Force f_atom = (*f2);
361 f[localID.
pid][localID.
index] += f_atom;
367 extVirial +=
outer(f_atom,x_atom);
370 DebugM(1,
"done with the loop\n");
374 g_i = gdef.
begin(); g_e = gdef.
end();
384 for ( ; g_i != g_e; ++g_i, ++gf_i ) {
387 for ( ; *g_i != -1; ++g_i ) {
389 localID = atomMap->
localID(*g_i);
391 lidx = localID.
index;
392 if ( lpid ==
notUsed || ! f[lpid] )
continue;
393 f_atom = accel * mass_soa[lpid][lidx];
396 CkPrintf(
"NAMD3-recv: group %d, Before Force (%8.6f, %8.6f, %8.6f) \n",
397 *g_i, force_soa_x[lpid][lidx], force_soa_y[lpid][lidx], force_soa_z[lpid][lidx]);
398 CkPrintf(
"NAMD3-recv: group %d, Added Force (%8.6f, %8.6f, %8.6f) \n", *g_i, f_atom.
x, f_atom.
y, f_atom.
z);
401 force_soa_x[lpid][lidx] += f_atom.
x;
402 force_soa_y[lpid][lidx] += f_atom.
y;
403 force_soa_z[lpid][lidx] += f_atom.
z;
404 #ifndef USE_GLOBALMASTER_VIRIAL_KERNEL 405 x_orig.
x = pos_soa_x[lpid][lidx];
406 x_orig.
y = pos_soa_y[lpid][lidx];
407 x_orig.
z = pos_soa_z[lpid][lidx];
408 trans.
i = transform_soa_i[lpid][lidx];
409 trans.
j = transform_soa_j[lpid][lidx];
410 trans.
k = transform_soa_k[lpid][lidx];
413 extVirial +=
outer(f_atom,x_atom);
418 for ( ; g_i != g_e; ++g_i, ++gf_i ) {
421 for ( ; *g_i != -1; ++g_i ) {
429 CkPrintf(
"NAMD2-recv: group %d, Before Force (%8.6f, %8.6f, %8.6f) \n",
431 CkPrintf(
"NAMD2-recv: group %d, Added Force (%8.6f, %8.6f, %8.6f) \n", *g_i, f_atom.
x, f_atom.
y, f_atom.
z);
434 f[localID.
pid][localID.
index] += f_atom;
440 extVirial +=
outer(f_atom,x_atom);
444 DebugM(1,
"done with the groups\n");
446 if (numActiveGridObjects > 0) {
447 applyGridObjectForces(msg, &extForce, &extVirial);
460 DebugM(3,
"Reconfiguring\n");
467 DebugM(3,
"Sending requested data right away\n");
472 groupTotalForce.
resize(numGroupsRequested);
473 for (
int i=0; i<numGroupsRequested; ++i ) groupTotalForce[i] = 0;
476 DebugM(3,
"setting forces\n");
478 Force **f = forcePtrs;
480 for (ap = ap.
begin(); ap != ap.
end(); ap++) {
484 (*ap).positionBox->close(&x);
485 (*ap).forceBox->close(&((*ap).r));
486 DebugM(1,
"closing boxes\n");
492 mass_soa[pId] = NULL;
493 pos_soa_x[pId] = NULL;
494 pos_soa_y[pId] = NULL;
495 pos_soa_z[pId] = NULL;
496 force_soa_x[pId] = NULL;
497 force_soa_y[pId] = NULL;
498 force_soa_z[pId] = NULL;
499 transform_soa_i[pId] = NULL;
500 transform_soa_j[pId] = NULL;
501 transform_soa_k[pId] = NULL;
502 DebugM(2,
"nulling ptrs\n");
505 DebugM(3,
"done setting forces\n");
508 #ifdef NODEGROUP_FORCE_REGISTER 516 DebugM(3,
"Done processing results\n");
521 DebugM(2,
"doWork thread " << CthGetToken(CthSelf())->serialNo <<
"\n");
529 DebugM(3,
"doWork for step " << step <<
"\n"<<
endi);
534 for (ap = ap.
begin(); ap != ap.
end(); ap++) {
535 CompAtom *x = (*ap).positionBox->open();
536 t[(*ap).patchID] = (*ap).p->getAtomList().begin();
541 mass_soa[pId] = (*ap).p->patchDataSOA.mass;
542 pos_soa_x[pId] = (*ap).p->patchDataSOA.pos_x;
543 pos_soa_y[pId] = (*ap).p->patchDataSOA.pos_y;
544 pos_soa_z[pId] = (*ap).p->patchDataSOA.pos_z;
545 transform_soa_i[pId] = (*ap).p->patchDataSOA.transform_i;
546 transform_soa_j[pId] = (*ap).p->patchDataSOA.transform_j;
547 transform_soa_k[pId] = (*ap).p->patchDataSOA.transform_k;
569 #ifdef NODEGROUP_FORCE_REGISTER 573 comm->stowSuspendULT();
576 comm->stowSuspendULT();
583 #endif // NODEGROUP_FORCE_REGISTER 593 DebugM(2,
"skipping step "<< step <<
"\n"<<
endi);
601 DebugM(2,
"done with doWork\n");
604 void ComputeGlobal::sendData()
632 for ( ; a != a_e; ++a ) {
633 localID = atomMap->
localID(*a);
635 lidx = localID.
index;
636 if ( lpid ==
notUsed || ! t[lpid] )
continue;
639 x_orig.
x = pos_soa_x[lpid][lidx];
640 x_orig.
y = pos_soa_y[lpid][lidx];
641 x_orig.
z = pos_soa_z[lpid][lidx];
642 trans.
i = transform_soa_i[lpid][lidx];
643 trans.
j = transform_soa_j[lpid][lidx];
644 trans.
k = transform_soa_k[lpid][lidx];
649 for ( ; a != a_e; ++a ) {
665 g_i = gdef.
begin(); g_e = gdef.
end();
672 for ( ; g_i != g_e; ++g_i ) {
675 for ( ; *g_i != -1; ++g_i ) {
676 localID = atomMap->
localID(*g_i);
678 lidx = localID.
index;
679 if ( lpid ==
notUsed || ! t[lpid] )
continue;
681 x_orig.
x = pos_soa_x[lpid][lidx];
682 x_orig.
y = pos_soa_y[lpid][lidx];
683 x_orig.
z = pos_soa_z[lpid][lidx];
684 trans.
i = transform_soa_i[lpid][lidx];
685 trans.
j = transform_soa_j[lpid][lidx];
686 trans.
k = transform_soa_k[lpid][lidx];
688 mass += mass_soa[lpid][lidx];
691 printf(
"NAMD3-send: step %d atom %d, POS (%8.6f, %8.6f, %8.6f) \n",
patchList[0].p->flags.step, *g_i, x_orig.
x, x_orig.
y, x_orig.
z);
698 DebugM(1,
"Adding center of mass "<<com<<
"\n");
705 for ( ; g_i != g_e; ++g_i ) {
708 for ( ; *g_i != -1; ++g_i ) {
719 printf(
"NAMD2-send: step %d atom %d, POS (%8.6f, %8.6f, %8.6f) \n",
patchList[0].p->flags.step, *g_i, x_orig.
x, x_orig.
y, x_orig.
z);
727 DebugM(1,
"Adding center of mass "<<com<<
"\n");
738 if (numActiveGridObjects > 0) {
740 computeGridObjects(msg);
748 if ( gfcount ) msg->
gtf.
swap(groupTotalForce);
752 DebugM(3,
"Sending data (" << msg->
p.
size() <<
" positions, " 754 <<
" grid objects) on client\n");
765 #ifdef NODEGROUP_FORCE_REGISTER 771 comm->stowSuspendULT();
772 comm->stowSuspendULT();
797 for ( ; ai != ae; ai++, gfi++) {
798 *gfi =
Vector(0.0, 0.0, 0.0);
808 DebugM(1,
"id = " << ai->
id <<
", scale = " << scale
809 <<
", charge = " << charge <<
", position = " << pos <<
"\n");
810 if (grid->compute_VdV(pos, V, dV)) {
815 *gfi = -charge * scale * dV;
816 gridObjValue += charge * scale * V;
817 DebugM(1,
"id = " << ai->
id <<
", force = " << *gfi <<
"\n");
819 DebugM(3,
"gridObjValue = " << gridObjValue <<
"\n" <<
endi);
829 NAMD_bug(
"No grids loaded in memory but ComputeGlobal has been requested to use them.");
837 size_t ig = 0, gridobjcount = 0;
841 for (ap = ap.begin(); ap != ap.end(); ap++) {
845 int const numAtoms = ap->p->getNumAtoms();
846 ForceList **gridForces = gridForcesPtrs[ap->p->getPatchID()];
849 for (ig = 0; ig < numGridObjects; ig++) {
851 DebugM(2,
"Processing grid index " << ig <<
"\n" <<
endi);
854 if (!gridObjActive[ig]) {
855 DebugM(2,
"Skipping grid index " << ig <<
"; it is handled by " 856 "ComputeGridForce\n" <<
endi);
860 ForceList *gridForcesGrid = gridForces[ig];
861 gridForcesGrid->
resize(numAtoms);
867 DebugM(2,
"computeGridObjects(): patch = " << ap->p->getPatchID()
868 <<
", grid index = " << ig <<
"\n" <<
endi);
874 computeGridForceGrid(ai, ae, gfi, ap->p->lattice, ig, grid, gridobjvalue);
880 for (gridobjcount = 0; gridobjcount < numActiveGridObjects; gridobjcount++) {
885 DebugM(2,
"computeGridObjects done\n");
895 NAMD_bug(
"ComputeGlobal received a different number of grid forces than active grids.");
901 Force &extForce = *extForce_in;
902 Tensor &extVirial = *extVirial_in;
906 gridObjForces.
resize(numGridObjects);
907 gridObjForces.
setall(0.0);
911 for (ig = 0; gridobjforce_i != gridobjforce_e ;
912 gridobjforce_i++, ig++) {
913 if (!gridObjActive[ig])
continue;
914 gridObjForces[ig] = *gridobjforce_i;
919 for (ap = ap.begin(); ap != ap.end(); ap++) {
921 ForceList **gridForces = gridForcesPtrs[ap->p->getPatchID()];
923 for (ig = 0; ig < numGridObjects; ig++) {
925 if (!gridObjActive[ig])
continue;
927 DebugM(2,
"gof = " << gridObjForces[ig] <<
"\n" <<
endi);
929 ForceList *gridForcesGrid = gridForces[ig];
936 for ( ; ai != ae; ai++, gfi++) {
938 *gfi =
Vector(0.0, 0.0, 0.0);
943 Vector const gridforceatom(-1.0 * (*gfi) * gridObjForces[ig]);
945 <<
", pid = " << localID.
pid 946 <<
", index = " << localID.
index 947 <<
", force = " << gridforceatom <<
"\n" <<
endi);
948 f[localID.
index] += gridforceatom;
949 extForce += gridforceatom;
953 extVirial +=
outer(gridforceatom, x_virial);
970 if ( ! forceSendEnabled )
NAMD_bug(
"ComputeGlobal::saveTotalForces called unexpectedly");
971 if ( ! forceSendActive )
return;
979 for (
int i=0; i<num; ++i) {
980 int index = atoms[i].id;
981 if (index < endRequested && isRequested[index] & 1) {
983 totalForce.
add(af[i]);
997 double *f1_soa_x = homePatch->patchDataSOA.
f_normal_x;
998 double *f1_soa_y = homePatch->patchDataSOA.
f_normal_y;
999 double *f1_soa_z = homePatch->patchDataSOA.
f_normal_z;
1008 double f_sum_x, f_sum_y, f_sum_z;
1011 for (
int i=0; i<num; ++i) {
1012 int index = atoms[i].id;
1015 CkPrintf(
"ForceSaved: atom %d, ForceN (%8.6f, %8.6f, %8.6f) \n", index, f1_soa_x[i], f1_soa_y[i], f1_soa_z[i]);
1016 CkPrintf(
" atom %d, ForceNB (%8.6f, %8.6f, %8.6f) \n", index, f2_soa_x[i], f2_soa_y[i], f2_soa_z[i]);
1017 CkPrintf(
" atom %d, ForceSL (%8.6f, %8.6f, %8.6f) \n", index, f3_soa_x[i], f3_soa_y[i], f3_soa_z[i]);
1019 CkPrintf(
"ForceSaved: atom %d, ForceN (%8.6f, %8.6f, %8.6f) \n", index, f1[i].x, f1[i].y, f1[i].z);
1020 CkPrintf(
" atom %d, ForceNB (%8.6f, %8.6f, %8.6f) \n", index, f2[i].x, f2[i].y, f2[i].z);
1027 printf(
"PE, PId (%d, %d) Stop saving at step: %d ####################################################\n",
1030 if ( ! forceSendActive )
return;
1031 for (
int i=0; i<num; ++i) {
1032 int index = atoms[i].id;
1034 if (index < endRequested && (reqflag = isRequested[index])) {
1036 f_sum_x = f1_soa_x[i] + f2_soa_x[i];
1037 f_sum_y = f1_soa_y[i] + f2_soa_y[i];
1038 f_sum_z = f1_soa_z[i] + f2_soa_z[i];
1040 f_sum_x += f3_soa_x[i];
1041 f_sum_y += f3_soa_y[i];
1042 f_sum_z += f3_soa_z[i];
1048 f_sum = f1[i]+f2[i];
1053 if ( fixedAtomsOn && atoms[i].atomFixed )
1056 if ( reqflag & 1 ) {
1058 totalForce.
add(f_sum);
1060 if ( reqflag & 2 ) {
1063 if ( gpi == gpend || gpi->
first != index )
1064 NAMD_bug(
"ComputeGlobal::saveTotalForces gpair corrupted.");
1067 groupTotalForce[gpi->
second] += f_sum;
1068 }
while ( ++gpi != gpend && gpi->
first == index );
#define NAMD_EVENT_STOP(eon, id)
GridforceGrid * get_gridfrc_grid(int gridnum) const
void recvResults(ComputeGlobalResultsMsg *)
NAMD_HOST_DEVICE Position reverse_transform(Position data, const Transform &t) const
Bool globalMasterScaleByFrequency
static PatchMap * Object()
BigRealList gridobjvalue
Partial values of the GridForce objects from this message.
NAMD_HOST_DEVICE Tensor outer(const Vector &v1, const Vector &v2)
virtual void submit(void)=0
#define ADD_TENSOR_OBJECT(R, RL, D)
SimParameters * simParameters
ComputeHomePatchList patchList
Bool CUDASOAintegrateMode
std::ostream & endi(std::ostream &s)
SubmitReduction * willSubmit(int setID, int size=-1)
void saveTotalForces(HomePatch *)
ResizeArrayIter< T > begin(void) const
static ReductionMgr * Object(void)
ComputeGlobal(ComputeID, ComputeMgr *)
int add(const Elem &elem)
Molecule stores the structural information for the system.
void setall(const Elem &elem)
ResizeArray< Force > ForceList
int numPatches(void) const
#define NAMD_EVENT_START(eon, id)
const Elem * const_iterator
void NAMD_bug(const char *err_msg)
LocalID localID(AtomID id)
ResizeArray< Lattice > lat
IntList gridobjindex
Indices of the GridForce objects contained in this message.
void get_gridfrc_params(Real &k, Charge &q, int atomnum, int gridnum) const
void enableComputeGlobalResults()
void recvComputeGlobalResults(ComputeGlobalResultsMsg *)
static AtomMap * Object()
#define ADD_VECTOR_OBJECT(R, RL, D)
int count
Numer of atoms processed for this message.
Bool is_atom_gridforced(int atomnum, int gridnum) const
void sendComputeGlobalData(ComputeGlobalDataMsg *)
ForceList f[Results::maxNumForces]
int patchcount
Number of patches processed for this message.
void swap(ResizeArray< Elem > &ra)
ResizeArrayIter< T > end(void) const
int globalMasterFrequency