From: 卢禹锟 (luyukun_at_itp.ac.cn)
Date: Sun Jun 15 2014 - 04:57:43 CDT
Dear namd developer and user,
Recently I try to understand the namd source code about the nonbonded force calculation, especially how it send the force to the integrator.
To my best knowledge, the nonbonded force types were first announced in the NamdTypes.h defined in the Results class in the way:
enum { normal=0, nbond=1, slow=2, amdf=3,maxNumForces=4 }
after some Initialization then it will be first called by the ComputeNonbondedPair.C and ComputeNonbondedSelf.C in the way :
params.ff[0] = r[a]->f[Results::nbond];
params.ff[1] = r[b]->f[Results::nbond];
.....
params.fullf[0] = r[a]->f[Results::slow];
params.fullf[1] = r[b]->f[Results::slow];
and followed by continuously adding the vdw and elect force pairs calculated from different patches in
ComputeNonbondedBase.h and ComputeNonbondedBase2.h (I will show it in a equivalent but simpler code) :
Force *f_0 = params->ff[0];
Force *f_1 = params->ff[1];
.....
Force *fullf_0 = params->fullf[0];
Force *fullf_1 = params->fullf[1];
.....
f_0[i].x += force_r * p_ij_x;
f_1[j].x -= force_r * p_ij_x;
......
fullf_0[i].x += fullforce_r * p_ij_x;
fullf_1[j].x -= fullforce_r * p_ij_x;
.....
When I only allow prescribed atom pairs to contribute to f[Results::nbond] and f[Results::slow], and ignore the pme contribution. Hence I modify ComputeNonbondedBase.h and ComputeNonbondedBase2.h to:
if (isSelPairs==1) f_0[i].x += force_r * p_ij_x; // isSelPairs=1 if the pairs is allowed
if (isSelPairs==1) f_1[j].x -= force_r * p_ij_x;
if (isSelPairs==1) fullf_0[i].x += fullforce_r * p_ij_x;
if (isSelPairs==1) fullf_1[j].x -= fullforce_r * p_ij_x;
and clear the pme force in ComputePme.C by:
Force *f = r->f[Results::slow];
for (int i = 0; i < numAtoms; ++i) {
f[i].x += results_ptr->x*0.; // multiply 0 to clear pme force
f[i].y += results_ptr->y*0.;
f[i].z += results_ptr->z*0.;
++results_ptr;
}
Then for the two atom pairs, the sum of f[Results::nbond] and f[Results::slow] (No matter how the nonbondedFreq and fullElectFrequency you set in the Tcl script) should be the short distance nonbonded force between them. If we print the pair force message in the integrate() after the force calculation function runComputeObjects() (all functions are defined in sequencer.C), the log file show the following results (Fx Fy Fz is the force in xyz direction):
SSSS: replicaWeight: 1 rescalingFactor: 1 DE: 31733.6 Fx: 0 Fy: 0 Fz: 0
SSSS: replicaWeight: 1 rescalingFactor: 1 DE: 31733.6 Fx: 1.3565 Fy: -0.592703 Fz: 1.17894
TIMING: 99 CPU: 2.36917, 0.022063/step Wall: 2.36917, 0.022063/step, 6.12862e-06 hours remaining, 3.646919 MB of memory in use.
ENERGY: 99 7.4593 21.9494 33.5562 0.7441 -12503.0952 1526.5896 0.0000 0.0000 2700.4726 -8212.3241 450.4275 -10912.7967 -8212.3241 450.4275 4896.1036 1915.9119 33600.0000 4896.1036 1915.9119
LDB: ============= START OF LOAD BALANCING ============== 2.42457
WRITING COORDINATES TO DCD FILE chargeala12.dcd AT STEP 99
LDB: Largest compute 3 load 0.090234 is 19.3% of average load 0.466756
LDB: Average compute 0.030192 is 6.5% of average load 0.466756
LDB: TIME 2.42511 LOAD: AVG 0.466756 MAX 0.908014 PROXIES: TOTAL 4 MAXPE 1 MAXPATCH 2 None MEM: 3.66975 MB
LDB: TIME 2.42518 LOAD: AVG 0.466756 MAX 0.468264 PROXIES: TOTAL 4 MAXPE 1 MAXPATCH 2 TorusLB MEM: 3.67107 MB
LDB: ============== END OF LOAD BALANCING =============== 2.42583
Info: useSync: 1 useProxySync: 0
LDB: =============== DONE WITH MIGRATION ================ 2.42715
SSSS: replicaWeight: 1 rescalingFactor: 1 DE: 31733.6 Fx: 0 Fy: 0 Fz: 0
SSSS: replicaWeight: 1 rescalingFactor: 1 DE: 31733.6 Fx: 0 Fy: 0 Fz: 0
The line started with "SSSS:" is the pairs force output message ( here namd will use two patches, so there are two output messages and only one patch contain the atom pairs I want, then the other one will have a zero force). After the load balancing, the force change from a correct value Fx: 1.3565 Fy: -0.592703 Fz: 1.17894 to Fx: 0 Fy: 0 Fz: 0.
I have try different pairs, wait for long simulation time to see if the force will change to nonzero value after several load balancing, but it show the same problem once after the first load balancing.
Could someone tell me why?
This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:22:28 CST