diff --git a/src/KSPACE/pair_lj_charmm_coul_long.cpp b/src/KSPACE/pair_lj_charmm_coul_long.cpp index 2c709a3ae0..c9ceaa92df 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_long.cpp @@ -172,9 +172,9 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -206,7 +206,7 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; evdwl *= switch1; } evdwl *= factor_lj; @@ -247,13 +247,6 @@ void PairLJCharmmCoulLong::compute_inner() numneigh = listinner->numneigh; firstneigh = listinner->firstneigh; - double cut_out_on = cut_respa[0]; - double cut_out_off = cut_respa[1]; - - double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { @@ -289,7 +282,7 @@ void PairLJCharmmCoulLong::compute_inner() fpair = (forcecoul + factor_lj*forcelj) * r2inv; if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + rsw = (sqrt(rsq) - cut_out_on)*cut_out_diff_inv; fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); } @@ -332,18 +325,6 @@ void PairLJCharmmCoulLong::compute_middle() numneigh = listmiddle->numneigh; firstneigh = listmiddle->firstneigh; - double cut_in_off = cut_respa[0]; - double cut_in_on = cut_respa[1]; - double cut_out_on = cut_respa[2]; - double cut_out_off = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { @@ -378,20 +359,20 @@ void PairLJCharmmCoulLong::compute_middle() forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } fpair = (forcecoul + factor_lj*forcelj) * r2inv; if (rsq < cut_in_on_sq) { - rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; + rsw = (sqrt(rsq) - cut_in_off)*cut_in_diff_inv; fpair *= rsw*rsw*(3.0 - 2.0*rsw); } if (rsq > cut_out_on_sq) { - rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; + rsw = (sqrt(rsq) - cut_out_on)*cut_out_diff_inv; fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); } @@ -441,13 +422,6 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) numneigh = listouter->numneigh; firstneigh = listouter->firstneigh; - double cut_in_off = cut_respa[2]; - double cut_in_on = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { @@ -518,9 +492,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -562,7 +536,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; evdwl *= switch1; } evdwl *= factor_lj; @@ -590,9 +564,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -601,9 +575,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -759,14 +733,27 @@ void PairLJCharmmCoulLong::init_style() cut_coulsq = cut_coul * cut_coul; cut_bothsq = MAX(cut_ljsq,cut_coulsq); - denom_lj = (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) * - (cut_ljsq-cut_lj_innersq); + denom_lj_inv = 1.0 / ( (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) * + (cut_ljsq-cut_lj_innersq) ); // set & error check interior rRESPA cutoffs if (strstr(update->integrate_style,"respa") && ((Respa *) update->integrate)->level_inner >= 0) { cut_respa = ((Respa *) update->integrate)->cutoff; + cut_in_off = cut_respa[0]; + cut_in_on = cut_respa[1]; + cut_out_on = cut_respa[2]; + cut_out_off = cut_respa[3]; + + cut_in_diff = cut_in_on - cut_in_off; + cut_out_diff = cut_out_off - cut_out_on; + cut_in_diff_inv = 1.0 / (cut_in_diff); + cut_out_diff_inv = 1.0 / (cut_out_diff); + cut_in_off_sq = cut_in_off*cut_in_off; + cut_in_on_sq = cut_in_on*cut_in_on; + cut_out_on_sq = cut_out_on*cut_out_on; + cut_out_off_sq = cut_out_off*cut_out_off; if (MIN(cut_lj,cut_coul) < cut_respa[3]) error->all(FLERR,"Pair cutoff < Respa interior cutoff"); if (cut_lj_inner < cut_respa[1]) @@ -992,9 +979,9 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype, forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -1017,7 +1004,7 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype, philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; philj *= switch1; } eng += factor_lj*philj; diff --git a/src/KSPACE/pair_lj_charmm_coul_long.h b/src/KSPACE/pair_lj_charmm_coul_long.h index 8e2c049088..1544f3bc14 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.h +++ b/src/KSPACE/pair_lj_charmm_coul_long.h @@ -54,7 +54,11 @@ class PairLJCharmmCoulLong : public Pair { double cut_lj_innersq,cut_ljsq; double cut_coul,cut_coulsq; double cut_bothsq; - double denom_lj; + double cut_in_off, cut_in_on, cut_out_off, cut_out_on; + double cut_in_diff, cut_out_diff; + double cut_in_diff_inv, cut_out_diff_inv; + double cut_in_off_sq, cut_in_on_sq, cut_out_off_sq, cut_out_on_sq; + double denom_lj, denom_lj_inv; double **epsilon,**sigma,**eps14,**sigma14; double **lj1,**lj2,**lj3,**lj4,**offset; double **lj14_1,**lj14_2,**lj14_3,**lj14_4; diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.cpp b/src/KSPACE/pair_lj_charmm_coul_msm.cpp index e99071f4ef..76c9ef0cc7 100644 --- a/src/KSPACE/pair_lj_charmm_coul_msm.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_msm.cpp @@ -161,9 +161,9 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -221,7 +221,7 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; evdwl *= switch1; } evdwl *= factor_lj; @@ -358,9 +358,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -402,7 +402,7 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; evdwl *= switch1; } evdwl *= factor_lj; @@ -430,9 +430,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -441,9 +441,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -499,9 +499,9 @@ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype, forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; switch2 = 12.0*rsq * (cut_ljsq-rsq) * - (rsq-cut_lj_innersq) / denom_lj; + (rsq-cut_lj_innersq) * denom_lj_inv; philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]); forcelj = forcelj*switch1 + philj*switch2; } @@ -524,7 +524,7 @@ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype, philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * - (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj; + (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv; philj *= switch1; } eng += factor_lj*philj; diff --git a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp index 48e23d5833..0a4ddd7908 100755 --- a/src/USER-MISC/pair_lj_sf_dipole_sf.cpp +++ b/src/USER-MISC/pair_lj_sf_dipole_sf.cpp @@ -133,8 +133,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) if (rsq < cut_coulsq[itype][jtype]) { + rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; if (qtmp != 0.0 && q[j] != 0.0) { - pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]); + pre1 = qtmp*q[j]*rinv*(r2inv-rcutcoul2inv); forcecoulx += pre1*delx; forcecouly += pre1*dely; @@ -144,7 +145,6 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) if (mu[i][3] > 0.0 && mu[j][3] > 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; - rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2]; pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; @@ -156,7 +156,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) aforcecouly = pre1*dely; aforcecoulz = pre1*delz; - bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) + + bfac = 1.0 - 4.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv + 3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv; presf = 2.0 * r2inv * pidotr * pjdotr; bforcecoulx = bfac * (pjdotr*mu[i][0]+pidotr*mu[j][0]-presf*delx); @@ -187,10 +187,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz; - rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pre1 = 3.0 * q[j] * r5inv * pidotr * (1-rsq*rcutcoul2inv); pqfac = 1.0 - 3.0*rsq*rcutcoul2inv + - 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); + 2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv; pre2 = q[j] * r3inv * pqfac; forcecoulx += pre2*mu[i][0] - pre1*delx; @@ -205,10 +204,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz; - rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pre1 = 3.0 * qtmp * r5inv * pjdotr * (1-rsq*rcutcoul2inv); qpfac = 1.0 - 3.0*rsq*rcutcoul2inv + - 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); + 2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv; pre2 = qtmp * r3inv * qpfac; forcecoulx += pre1*delx - pre2*mu[j][0]; @@ -261,7 +259,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) if (eflag) { if (rsq < cut_coulsq[itype][jtype]) { - ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype])); + ecoul = (1.0-sqrt(rsq/cut_coulsq[itype][jtype])); ecoul *= ecoul; ecoul *= qtmp * q[j] * rinv; if (mu[i][3] > 0.0 && mu[j][3] > 0.0) diff --git a/src/USER-OMP/fix_nh_sphere_omp.cpp b/src/USER-OMP/fix_nh_sphere_omp.cpp index 71c5e4b985..5e7890a513 100644 --- a/src/USER-OMP/fix_nh_sphere_omp.cpp +++ b/src/USER-OMP/fix_nh_sphere_omp.cpp @@ -55,7 +55,7 @@ void FixNHSphereOMP::init() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (radius[i] == 0.0) - error->one(FLERR,"Fix nvt/sphere requires extended particles"); + error->one(FLERR,"Fix nvt/npt/nph/sphere/omp require extended particles"); FixNHOMP::init(); } @@ -97,6 +97,7 @@ void FixNHSphereOMP::nve_v() v[i].x += dtfm*f[i].x; v[i].y += dtfm*f[i].y; v[i].z += dtfm*f[i].z; + const double dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); omega[i].x += dtirotate*torque[i].x; omega[i].y += dtirotate*torque[i].y; diff --git a/src/USER-OMP/fix_nve_sphere_omp.cpp b/src/USER-OMP/fix_nve_sphere_omp.cpp index ee58c88cb8..2367e34878 100644 --- a/src/USER-OMP/fix_nve_sphere_omp.cpp +++ b/src/USER-OMP/fix_nve_sphere_omp.cpp @@ -21,13 +21,17 @@ #include "respa.h" #include "force.h" #include "error.h" +#include "math_vector.h" +#include "math_extra.h" using namespace LAMMPS_NS; using namespace FixConst; +using namespace MathExtra; #define INERTIA 0.4 // moment of inertia prefactor for sphere enum{NONE,DIPOLE}; +enum{NODLM,DLM}; /* ---------------------------------------------------------------------- */ @@ -77,21 +81,128 @@ void FixNVESphereOMP::initial_integrate(int vflag) if (extra == DIPOLE) { double * const * const mu = atom->mu; + if (dlm == NODLM) { #if defined(_OPENMP) #pragma omp parallel for private(i) default(none) #endif - for (i = 0; i < nlocal; i++) { - double g0,g1,g2,msq,scale; - if (mask[i] & groupbit) { - if (mu[i][3] > 0.0) { - g0 = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); - g1 = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); - g2 = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); - msq = g0*g0 + g1*g1 + g2*g2; - scale = mu[i][3]/sqrt(msq); - mu[i][0] = g0*scale; - mu[i][1] = g1*scale; - mu[i][2] = g2*scale; + for (i = 0; i < nlocal; i++) { + double g0,g1,g2,msq,scale; + if (mask[i] & groupbit) { + if (mu[i][3] > 0.0) { + g0 = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); + g1 = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); + g2 = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); + msq = g0*g0 + g1*g1 + g2*g2; + scale = mu[i][3]/sqrt(msq); + mu[i][0] = g0*scale; + mu[i][1] = g1*scale; + mu[i][2] = g2*scale; + } + } + } + } else { +#if defined(_OPENMP) +#pragma omp parallel for private(i) default(none) +#endif + // Integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme + for (i = 0; i < nlocal; i++) { + vector w, w_temp, a; + matrix Q, Q_temp, R; + + if (mask[i] & groupbit && mu[i][3] > 0.0) { + + // Construct Q from dipole: + // Q is the rotation matrix from space frame to body frame + // i.e. v_b = Q.v_s + + // Define mu to lie along the z axis in the body frame + // We take the unit dipole to avoid getting a scaling matrix + const double inv_len_mu = 1.0/mu[i][3]; + a[0] = mu[i][0]*inv_len_mu; + a[1] = mu[i][1]*inv_len_mu; + a[2] = mu[i][2]*inv_len_mu; + + // v = a x [0 0 1] - cross product of mu in space and body frames + // s = |v| + // c = a.[0 0 1] = a[2] + // vx = [ 0 -v[2] v[1] + // v[2] 0 -v[0] + // -v[1] v[0] 0 ] + // then + // Q = I + vx + vx^2 * (1-c)/s^2 + + const double s2 = a[0]*a[0] + a[1]*a[1]; + if (s2 != 0.0){ // i.e. the vectors are not parallel + const double scale = (1.0 - a[2])/s2; + + Q[0][0] = 1.0 - scale*a[0]*a[0]; Q[0][1] = -scale*a[0]*a[1]; Q[0][2] = -a[0]; + Q[1][0] = -scale*a[0]*a[1]; Q[1][1] = 1.0 - scale*a[1]*a[1]; Q[1][2] = -a[1]; + Q[2][0] = a[0]; Q[2][1] = a[1]; Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]); + } else { // if parallel then we just have I or -I + Q[0][0] = 1.0/a[2]; Q[0][1] = 0.0; Q[0][2] = 0.0; + Q[1][0] = 0.0; Q[1][1] = 1.0/a[2]; Q[1][2] = 0.0; + Q[2][0] = 0.0; Q[2][1] = 0.0; Q[2][2] = 1.0/a[2]; + } + + // Local copy of this particle's angular velocity (in space frame) + w[0] = omega[i][0]; w[1] = omega[i][1]; w[2] = omega[i][2]; + + // Transform omega into body frame: w_temp= Q.w + matvec(Q,w,w_temp); + + // Construct rotation R1 + BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]); + + // Apply R1 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R1 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Construct rotation R2 + BuildRyMatrix(R, dtf/force->ftm2v*w[1]); + + // Apply R2 to w: w_temp = R.w + matvec(R,w,w_temp); + + // Apply R2 to Q: Q = R^T.Q_temp + transpose_times3(R,Q_temp,Q); + + // Construct rotation R3 + BuildRzMatrix(R, 2.0*dtf/force->ftm2v*w_temp[2]); + + // Apply R3 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R3 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Construct rotation R4 + BuildRyMatrix(R, dtf/force->ftm2v*w[1]); + + // Apply R4 to w: w_temp = R.w + matvec(R,w,w_temp); + + // Apply R4 to Q: Q = R^T.Q_temp + transpose_times3(R,Q_temp,Q); + + // Construct rotation R5 + BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]); + + // Apply R5 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R5 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Transform w back into space frame w_temp = Q^T.w + transpose_matvec(Q_temp,w,w_temp); + omega[i][0] = w_temp[0]; omega[i][1] = w_temp[1]; omega[i][2] = w_temp[2]; + + // Set dipole according to updated Q: mu = Q^T.[0 0 1] * |mu| + mu[i][0] = Q_temp[2][0] * mu[i][3]; + mu[i][1] = Q_temp[2][1] * mu[i][3]; + mu[i][2] = Q_temp[2][2] * mu[i][3]; } } } diff --git a/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp b/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp index 67afbbeefe..1eadfe3d1b 100755 --- a/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp +++ b/src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp @@ -123,8 +123,8 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - factor_coul = special_coul[sbmask(j)]; factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; j &= NEIGHMASK; delx = xtmp - x[j].x; @@ -137,8 +137,6 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) r2inv = 1.0/rsq; rinv = sqrt(r2inv); - // atom can have both a charge and dipole - // i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole // atom can have both a charge and dipole // i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole @@ -148,8 +146,9 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) if (rsq < cut_coulsq[itype][jtype]) { + rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; if (qtmp != 0.0 && q[j] != 0.0) { - pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]); + pre1 = qtmp*q[j]*rinv*(r2inv-rcutcoul2inv); forcecoulx += pre1*delx; forcecouly += pre1*dely; @@ -159,7 +158,6 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) if (mu[i].w > 0.0 && mu[j].w > 0.0) { r3inv = r2inv*rinv; r5inv = r3inv*r2inv; - rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pdotp = mu[i].x*mu[j].x + mu[i].y*mu[j].y + mu[i].z*mu[j].z; pidotr = mu[i].x*delx + mu[i].y*dely + mu[i].z*delz; @@ -171,7 +169,7 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) aforcecouly = pre1*dely; aforcecoulz = pre1*delz; - bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) + + bfac = 1.0 - 4.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv + 3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv; presf = 2.0 * r2inv * pidotr * pjdotr; bforcecoulx = bfac * (pjdotr*mu[i].x+pidotr*mu[j].x-presf*delx); @@ -202,10 +200,9 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pidotr = mu[i].x*delx + mu[i].y*dely + mu[i].z*delz; - rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pre1 = 3.0 * q[j] * r5inv * pidotr * (1-rsq*rcutcoul2inv); pqfac = 1.0 - 3.0*rsq*rcutcoul2inv + - 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); + 2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv; pre2 = q[j] * r3inv * pqfac; forcecoulx += pre2*mu[i].x - pre1*delx; @@ -220,10 +217,9 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr) r3inv = r2inv*rinv; r5inv = r3inv*r2inv; pjdotr = mu[j].x*delx + mu[j].y*dely + mu[j].z*delz; - rcutcoul2inv=1.0/cut_coulsq[itype][jtype]; pre1 = 3.0 * qtmp * r5inv * pjdotr * (1-rsq*rcutcoul2inv); qpfac = 1.0 - 3.0*rsq*rcutcoul2inv + - 2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv); + 2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv; pre2 = qtmp * r3inv * qpfac; forcecoulx += pre1*delx - pre2*mu[j].x; diff --git a/src/USER-QMMM/fix_qmmm.cpp b/src/USER-QMMM/fix_qmmm.cpp index c724d7b9a1..50adf8d224 100644 --- a/src/USER-QMMM/fix_qmmm.cpp +++ b/src/USER-QMMM/fix_qmmm.cpp @@ -34,7 +34,8 @@ // have to match with those from the QM code enum {QMMM_TAG_OTHER=0, QMMM_TAG_SIZE=1, QMMM_TAG_COORD=2, QMMM_TAG_FORCE=3, QMMM_TAG_FORCE2=4, QMMM_TAG_CELL=5, - QMMM_TAG_RADII=6, QMMM_TAG_CHARGE=7}; + QMMM_TAG_RADII=6, QMMM_TAG_CHARGE=7, QMMM_TAG_TYPE=8, + QMMM_TAG_MASS=9}; using namespace LAMMPS_NS; using namespace FixConst; @@ -406,8 +407,10 @@ void FixQMMM::exchange_positions() double *charge = atom->q; const int * const mask = atom->mask; const int * const type = atom->type; + const double * const mass = atom->mass; const tagint * const tag = atom->tag; const int nlocal = atom->nlocal; + const int ntypes = atom->ntypes; const int natoms = (int) atom->natoms; if (atom->natoms > MAXSMALLINT) @@ -422,12 +425,8 @@ void FixQMMM::exchange_positions() - create hash table that connects master list index with tag - collect necessary data and get master list index via hash table */ - if (comm->me == 0) { // AK: make aradii a per-atom property managed by - // fix qmmm so it migrates with the local atoms. - double *aradii = (double *) calloc(natoms, sizeof(double)); - ec_fill_radii( aradii, natoms ); - - // This array will pack all details about the cell + if (comm->me == 0) { + // Pack various cell dimension properties into one chunk. double celldata[9]; celldata[0] = domain->boxlo[0]; celldata[1] = domain->boxlo[1]; @@ -440,14 +439,13 @@ void FixQMMM::exchange_positions() celldata[8] = domain->yz; if (qmmm_role == QMMM_ROLE_MASTER) { - int isend_buf[3]; + int isend_buf[4]; isend_buf[0] = natoms; isend_buf[1] = num_qm; - isend_buf[2] =num_mm; - MPI_Send( isend_buf, 3, MPI_INTEGER,1, QMMM_TAG_SIZE, qm_comm ); - MPI_Send( celldata, 9, MPI_DOUBLE, 1, QMMM_TAG_CELL, qm_comm ); - MPI_Send( aradii, natoms, MPI_DOUBLE, 1, QMMM_TAG_RADII, qm_comm ); - free(aradii); + isend_buf[2] = num_mm; + isend_buf[3] = ntypes; + MPI_Send(isend_buf, 4, MPI_INTEGER,1, QMMM_TAG_SIZE, qm_comm); + MPI_Send(celldata, 9, MPI_DOUBLE, 1, QMMM_TAG_CELL, qm_comm); } if (verbose > 0) { if (screen) fputs("QMMM: exchange positions\n",screen); @@ -497,6 +495,8 @@ void FixQMMM::exchange_positions() MPI_Send(mm_charge_all, natoms, MPI_DOUBLE, 1, QMMM_TAG_COORD, qm_comm); MPI_Send(mm_coord_all, 3*natoms, MPI_DOUBLE, 1, QMMM_TAG_COORD, qm_comm); MPI_Send(mm_mask_all, natoms, MPI_INT, 1, QMMM_TAG_COORD, qm_comm); + MPI_Send(type, natoms, MPI_INT, 1, QMMM_TAG_TYPE, qm_comm); + MPI_Send(mass, ntypes+1, MPI_DOUBLE, 1, QMMM_TAG_MASS, qm_comm); /* to MM slave code */ MPI_Send(qm_coord, 3*num_qm, MPI_DOUBLE, 1, QMMM_TAG_COORD, mm_comm); @@ -838,247 +838,6 @@ double FixQMMM::memory_usage(void) } -/* ---------------------------------------------------------------------- */ -/* Manage the atomic number */ - -#define QMMM_ISOTOPES 351 - -static const int FixQMMM_Z[QMMM_ISOTOPES] = { - 1, 1, 1, 2, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 8, 9, - 10, 10, 10, 11, 12, 12, 12, 13, 14, 14, 14, 15, 16, 16, 16, - 16, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20, 20, 20, - 21, 22, 22, 22, 22, 22, 23, 23, 24, 24, 24, 24, 25, 26, 26, - 26, 26, 27, 28, 28, 28, 28, 28, 29, 29, 30, 30, 30, 30, 30, - 31, 31, 32, 32, 32, 32, 32, 33, 34, 34, 34, 34, 34, 34, 35, - 35, 36, 36, 36, 36, 36, 36, 37, 37, 38, 38, 38, 38, 39, 40, - 40, 40, 40, 40, 41, 42, 42, 42, 42, 42, 42, 42, 43, 43, 43, - 44, 44, 44, 44, 44, 44, 44, 45, 46, 46, 46, 46, 46, 46, 47, - 47, 48, 48, 48, 48, 48, 48, 48, 48, 49, 49, 50, 50, 50, 50, - 50, 50, 50, 50, 50, 50, 51, 51, 52, 52, 52, 52, 52, 52, 52, - 52, 53, 54, 54, 54, 54, 54, 54, 54, 54, 54, 55, 56, 56, 56, - 56, 56, 56, 56, 57, 57, 58, 58, 58, 58, 59, 60, 60, 60, 60, - 60, 60, 60, 61, 61, 62, 62, 62, 62, 62, 62, 62, 63, 63, 64, - 64, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 66, 66, 66, 67, - 68, 68, 68, 68, 68, 68, 69, 70, 70, 70, 70, 70, 70, 70, 71, - 71, 72, 72, 72, 72, 72, 72, 73, 73, 74, 74, 74, 74, 74, 75, - 75, 76, 76, 76, 76, 76, 76, 76, 77, 77, 78, 78, 78, 78, 78, - 78, 79, 80, 80, 80, 80, 80, 80, 80, 81, 81, 82, 82, 82, 82, - 83, 84, 84, 85, 85, 86, 86, 86, 87, 88, 88, 88, 88, 89, 90, - 90, 91, 92, 92, 92, 92, 92, 93, 93, 94, 94, 94, 94, 94, 94, - 95, 95, 96, 96, 96, 96, 96, 96, 97, 97, 98, 98, 98, 98, 99, - 100, 101, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, - 111, 112, 113, 114, 115}; - -static const char FixQMMM_EL[QMMM_ISOTOPES][4] = { - "H", "H", "H", "He", "He", "Li", "Li", "Be", "B", "B", - "C", "C", "C", "N", "N", "O", "O", "O", "F", "Ne", "Ne", - "Ne", "Na", "Mg", "Mg", "Mg", "Al", "Si", "Si", "Si", - "P", "S", "S", "S", "S", "Cl", "Cl", "Ar", "Ar", "Ar", - "K", "K", "K", "Ca", "Ca", "Ca", "Ca", "Ca", "Ca", "Sc", - "Ti", "Ti", "Ti", "Ti", "Ti", "V", "V", "Cr", "Cr", "Cr", - "Cr", "Mn", "Fe", "Fe", "Fe", "Fe", "Co", "Ni", "Ni", - "Ni", "Ni", "Ni", "Cu", "Cu", "Zn", "Zn", "Zn", "Zn", "Zn", - "Ga", "Ga", "Ge", "Ge", "Ge", "Ge", "Ge", "As", "Se", - "Se", "Se", "Se", "Se", "Se", "Br", "Br", "Kr", "Kr", - "Kr", "Kr", "Kr", "Kr", "Rb", "Rb", "Sr", "Sr", "Sr", - "Sr", "Y", "Zr", "Zr", "Zr", "Zr", "Zr", "Nb", "Mo", "Mo", - "Mo", "Mo", "Mo", "Mo", "Mo", "Tc", "Tc", "Tc", "Ru", "Ru", - "Ru", "Ru", "Ru", "Ru", "Ru", "Rh", "Pd", "Pd", "Pd", "Pd", - "Pd", "Pd", "Ag", "Ag", "Cd", "Cd", "Cd", "Cd", "Cd", "Cd", - "Cd", "Cd", "In", "In", "Sn", "Sn", "Sn", "Sn", "Sn", "Sn", - "Sn", "Sn", "Sn", "Sn", "Sb", "Sb", "Te", "Te", "Te", "Te", - "Te", "Te", "Te", "Te", "I", "Xe", "Xe", "Xe", "Xe", "Xe", - "Xe", "Xe", "Xe", "Xe", "Cs", "Ba", "Ba", "Ba", "Ba", "Ba", - "Ba", "Ba", "La", "La", "Ce", "Ce", "Ce", "Ce", "Pr", "Nd", - "Nd", "Nd", "Nd", "Nd", "Nd", "Nd", "Pm", "Pm", "Sm", "Sm", - "Sm", "Sm", "Sm", "Sm", "Sm", "Eu", "Eu", "Gd", "Gd", "Gd", - "Gd", "Gd", "Gd", "Gd", "Tb", "Dy", "Dy", "Dy", "Dy", "Dy", - "Dy", "Dy", "Ho", "Er", "Er", "Er", "Er", "Er", "Er", "Tm", - "Yb", "Yb", "Yb", "Yb", "Yb", "Yb", "Yb", "Lu", "Lu", "Hf", - "Hf", "Hf", "Hf", "Hf", "Hf", "Ta", "Ta", "W", "W", "W", "W", - "W", "Re", "Re", "Os", "Os", "Os", "Os", "Os", "Os", "Os", - "Ir", "Ir", "Pt", "Pt", "Pt", "Pt", "Pt", "Pt", "Au", "Hg", - "Hg", "Hg", "Hg", "Hg", "Hg", "Hg", "Tl", "Tl", "Pb", "Pb", - "Pb", "Pb", "Bi", "Po", "Po", "At", "At", "Rn", "Rn", "Rn", - "Fr", "Ra", "Ra", "Ra", "Ra", "Ac", "Th", "Th", "Pa", "U", - "U", "U", "U", "U", "Np", "Np", "Pu", "Pu", "Pu", "Pu", "Pu", - "Pu", "Am", "Am", "Cm", "Cm", "Cm", "Cm", "Cm", "Cm", "Bk", - "Bk", "Cf", "Cf", "Cf", "Cf", "Es", "Fm", "Md", "Md", "No", - "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn", - "Uut", "Uuq", "Uup"}; - -static const int FixQMMM_A[QMMM_ISOTOPES] = { - 1, 2, 3, 3, 4, 6, 7, 9, 10, 11, 12, 13, 14, 14, 15, 16, 17, 18, - 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, - 36, 35, 37, 36, 38, 40, 39, 40, 41, 40, 42, 43, 44, 46, 48, 45, - 46, 47, 48, 49, 50, 50, 51, 50, 52, 53, 54, 55, 54, 56, 57, 58, - 59, 58, 60, 61, 62, 64, 63, 65, 64, 66, 67, 68, 70, 69, 71, 70, - 72, 73, 74, 76, 75, 74, 76, 77, 78, 80, 82, 79, 81, 78, 80, 82, - 83, 84, 86, 85, 87, 84, 86, 87, 88, 89, 90, 91, 92, 94, 96, 93, - 92, 94, 95, 96, 97, 98, 100, 97, 98, 99, 96, 98, 99, 100, 101, - 102, 104, 103, 102, 104, 105, 106, 108, 110, 107, 109, 106, 108, - 110, 111, 112, 113, 114, 116, 113, 115, 112, 114, 115, 116, 117, - 118, 119, 120, 122, 124, 121, 123, 120, 122, 123, 124, 125, 126, - 128, 130, 127, 124, 126, 128, 129, 130, 131, 132, 134, 136, 133, - 130, 132, 134, 135, 136, 137, 138, 138, 139, 136, 138, 140, 142, - 141, 142, 143, 144, 145, 146, 148, 150, 145, 147, 144, 147, 148, - 149, 150, 152, 154, 151, 153, 152, 154, 155, 156, 157, 158, 160, - 159, 156, 158, 160, 161, 162, 163, 164, 165, 162, 164, 166, 167, - 168, 170, 169, 168, 170, 171, 172, 173, 174, 176, 175, 176, 174, - 176, 177, 178, 179, 180, 180, 181, 180, 182, 183, 184, 186, 185, - 187, 184, 186, 187, 188, 189, 190, 192, 191, 193, 190, 192, 194, - 195, 196, 198, 197, 196, 198, 199, 200, 201, 202, 204, 203, 205, - 204, 206, 207, 208, 209, 209, 210, 210, 211, 211, 220, 222, 223, - 223, 224, 226, 228, 227, 230, 232, 231, 233, 234, 235, 236, 238, - 236, 237, 238, 239, 240, 241, 242, 244, 241, 243, 243, 244, 245, - 246, 247, 248, 247, 249, 249, 250, 251, 252, 252, 257, 258, 260, - 259, 262, 265, 268, 271, 272, 270, 276, 281, 280, 285, 284, 289, - 288}; - -static const double FixQMMM_MASS[QMMM_ISOTOPES] = { - 1.00782503207, 2.0141017778, 3.0160492777, 3.0160293191, 4.00260325415, 6.015122795, 7.01600455, 9.0121822, - 10.0129370, 11.0093054, 12.0000000, 13.0033548378, 14.003241989, 14.0030740048, 15.0001088982, 15.99491461956, - 16.99913170, 17.9991610, 18.99840322, 19.9924401754, 20.99384668, 21.991385114, 22.9897692809, 23.985041700, - 24.98583692, 25.982592929, 26.98153863, 27.9769265325, 28.976494700, 29.97377017, 30.97376163, 31.97207100, - 32.97145876, 33.96786690, 35.96708076, 34.96885268, 36.96590259, 35.967545106, 37.9627324, 39.9623831225, - 38.96370668, 39.96399848, 40.96182576, 39.96259098, 41.95861801, 42.9587666, 43.9554818, 45.9536926, - 47.952534, 44.9559119, 45.9526316, 46.9517631, 47.9479463, 48.9478700, 49.9447912, 49.9471585, - 50.9439595, 49.9460442, 51.9405075, 52.9406494, 53.9388804, 54.9380451, 53.9396105, 55.9349375, - 56.9353940, 57.9332756, 58.9331950, 57.9353429, 59.9307864, 60.9310560, 61.9283451, 63.9279660, - 62.9295975, 64.9277895, 63.9291422, 65.9260334, 66.9271273, 67.9248442, 69.9253193, 68.9255736, - 70.9247013, 69.9242474, 71.9220758, 72.9234589, 73.9211778, 75.9214026, 74.9215965, 73.9224764, - 75.9192136, 76.9199140, 77.9173091, 79.9165213, 81.9166994, 78.9183371, 80.9162906, 77.9203648, - 79.9163790, 81.9134836, 82.914136, 83.911507, 85.91061073, 84.911789738, 86.909180527, 83.913425, - 85.9092602, 86.9088771, 87.9056121, 88.9058483, 89.9047044, 90.9056458, 91.9050408, 93.9063152, - 95.9082734, 92.9063781, 91.906811, 93.9050883, 94.9058421, 95.9046795, 96.9060215, 97.9054082, - 99.907477, 96.906365, 97.907216, 98.9062547, 95.907598, 97.905287, 98.9059393, 99.9042195, - 100.9055821, 101.9043493, 103.905433, 102.905504, 101.905609, 103.904036, 104.905085, 105.903486, - 107.903892, 109.905153, 106.905097, 108.904752, 105.906459, 107.904184, 109.9030021, 110.9041781, - 111.9027578, 112.9044017, 113.9033585, 115.904756, 112.904058, 114.903878, 111.904818, 113.902779, - 114.903342, 115.901741, 116.902952, 117.901603, 118.903308, 119.9021947, 121.9034390, 123.9052739, - 120.9038157, 122.9042140, 119.904020, 121.9030439, 122.9042700, 123.9028179, 124.9044307, 125.9033117, - 127.9044631, 129.9062244, 126.904473, 123.9058930, 125.904274, 127.9035313, 128.9047794, 129.9035080, - 130.9050824, 131.9041535, 133.9053945, 135.907219, 132.905451933, 129.9063208, 131.9050613, 133.9045084, - 134.9056886, 135.9045759, 136.9058274, 137.9052472, 137.907112, 138.9063533, 135.907172, 137.905991, - 139.9054387, 141.909244, 140.9076528, 141.9077233, 142.9098143, 143.9100873, 144.9125736, 145.9131169, - 147.916893, 149.920891, 144.912749, 146.9151385, 143.911999, 146.9148979, 147.9148227, 148.9171847, - 149.9172755, 151.9197324, 153.9222093, 150.9198502, 152.9212303, 151.9197910, 153.9208656, 154.9226220, - 155.9221227, 156.9239601, 157.9241039, 159.9270541, 158.9253468, 155.924283, 157.924409, 159.9251975, - 160.9269334, 161.9267984, 162.9287312, 163.9291748, 164.9303221, 161.928778, 163.929200, 165.9302931, - 166.9320482, 167.9323702, 169.9354643, 168.9342133, 167.933897, 169.9347618, 170.9363258, 171.9363815, - 172.9382108, 173.9388621, 175.9425717, 174.9407718, 175.9426863, 173.940046, 175.9414086, 176.9432207, - 177.9436988, 178.9458161, 179.9465500, 179.9474648, 180.9479958, 179.946704, 181.9482042, 182.9502230, - 183.9509312, 185.9543641, 184.9529550, 186.9557531, 183.9524891, 185.9538382, 186.9557505, 187.9558382, - 188.9581475, 189.9584470, 191.9614807, 190.9605940, 192.9629264, 189.959932, 191.9610380, 193.9626803, - 194.9647911, 195.9649515, 197.967893, 196.9665687, 195.965833, 197.9667690, 198.9682799, 199.9683260, - 200.9703023, 201.9706430, 203.9734939, 202.9723442, 204.9744275, 203.9730436, 205.9744653, 206.9758969, - 207.9766521, 208.9803987, 208.9824304, 209.9828737, 209.987148, 210.9874963, 210.990601, 220.0113940, - 222.0175777, 223.0197359, 223.0185022, 224.0202118, 226.0254098, 228.0310703, 227.0277521, 230.0331338, - 232.0380553, 231.0358840, 233.0396352, 234.0409521, 235.0439299, 236.0455680, 238.0507882, 236.046570, - 237.0481734, 238.0495599, 239.0521634, 240.0538135, 241.0568515, 242.0587426, 244.064204, 241.0568291, - 243.0613811, 243.0613891, 244.0627526, 245.0654912, 246.0672237, 247.070354, 248.072349, 247.070307, - 249.0749867, 249.0748535, 250.0764061, 251.079587, 252.081626, 252.082980, 257.095105, 258.098431, - 260.10365, 259.10103, 262.10963, 265.11670, 268.12545, 271.13347, 272.13803, 270.13465, - 276.15116, 281.16206, 280.16447, 285.17411, 284.17808, 289.18728, 288.19249}; - -/* - * This function matches the element which has the least absolute - * difference between the masses. delta is set to the difference - * between the mass provided and the best match. - * - * THIS FUNCTION RELIES ON THE CORRECT ORDER OF THE DATA TABLES - * IN ORDER TO SKIP ISOTOPES FROM THE MATCH - * - * */ -int match_element(double mass, int search_isotopes, double &delta) -{ - int i; - int bestcandidate; - double diff, mindiff; - int lastz; - - mindiff = 1e6; - bestcandidate = -1; - lastz = -1; - for(i=0; icreate(type2element,atom->ntypes+1,"qmmm:type2element"); - type2element[0] = 0; - - for (i=1; i <= atom->ntypes; ++i) { - double delta; - el = match_element(atom->mass[i], 0, delta); - sprintf(errmsg,"FixQMMM: type %2d (mass: %8g) matches %2s with:" - " Z = %-3d A = %-3d r_cov = %5.2f (error = %-8.2g -> %-.2g%%)\n", - i, atom->mass[i], FixQMMM_EL[el], FixQMMM_Z[el], FixQMMM_A[el], - ec_r_covalent[FixQMMM_Z[el]], delta, delta/atom->mass[i] * 100.0); - type2element[i] = FixQMMM_Z[el]; - if(screen) fprintf(screen, "%s", errmsg); - if(logfile) fprintf(logfile, "%s", errmsg); - } - - for (i=0; itype[i]]; - if (el < 0 || el > EC_ELEMENTS) { - sprintf(errmsg,"Unable to find element Z=%d in table of covalent radii", el); - error->one(FLERR,errmsg); - } - aradii[i] = ec_r_covalent[el]; - if (ec_r_covalent[el] < 0.0) { - sprintf(errmsg, "Covalent radius for atom of element Z=%d not availabe", el); - error->one(FLERR,errmsg); - } - } - memory->destroy(type2element); - return 0; -} - // Local Variables: // mode: c++ // compile-command: "make -j4 openmpi" diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 37a813366c..8c576a8eee 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -79,6 +79,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) etap_mass_flag = 0; flipflag = 1; dipole_flag = 0; + dlm_flag = 0; tcomputeflag = 0; pcomputeflag = 0; @@ -330,7 +331,10 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) } else if (strcmp(arg[iarg],"update") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1; - else error->all(FLERR,"Illegal fix nvt/npt/nph command"); + else if (strcmp(arg[iarg],"dipole/dlm") == 0) { + dipole_flag = 1; + dlm_flag = 1; + } else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); @@ -434,7 +438,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (!atom->mu_flag) error->all(FLERR,"Using update dipole flag requires atom attribute mu"); } - + if ((tstat_flag && t_period <= 0.0) || (p_flag[0] && p_period[0] <= 0.0) || (p_flag[1] && p_period[1] <= 0.0) || diff --git a/src/fix_nh.h b/src/fix_nh.h index b0add1060b..c65014724a 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -111,6 +111,7 @@ class FixNH : public Fix { int omega_mass_flag; // 1 if omega_mass updated, 0 if not. int etap_mass_flag; // 1 if etap_mass updated, 0 if not. int dipole_flag; // 1 if dipole is updated, 0 if not. + int dlm_flag; // 1 if using the DLM rotational integrator, 0 if not int scaleyz; // 1 if yz scaled with lz int scalexz; // 1 if xz scaled with lz @@ -219,6 +220,10 @@ Self-explanatory. E: Using update dipole flag requires atom attribute mu +Self-explanatory. + +E: The dlm flag must be used with update dipole + Self-explanatory. E: Fix nvt/npt/nph damping parameters must be > 0.0 diff --git a/src/fix_nh_sphere.cpp b/src/fix_nh_sphere.cpp index f978b67660..de5ca0b3aa 100644 --- a/src/fix_nh_sphere.cpp +++ b/src/fix_nh_sphere.cpp @@ -21,9 +21,13 @@ #include "atom_vec.h" #include "group.h" #include "error.h" +#include "force.h" +#include "math_vector.h" +#include "math_extra.h" using namespace LAMMPS_NS; using namespace FixConst; +using namespace MathExtra; #define INERTIA 0.4 // moment of inertia prefactor for sphere @@ -50,7 +54,7 @@ void FixNHSphere::init() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (radius[i] == 0.0) - error->one(FLERR,"Fix nvt/sphere requires extended particles"); + error->one(FLERR,"Fix nvt/npt/nph/sphere require extended particles"); FixNH::init(); } @@ -102,28 +106,133 @@ void FixNHSphere::nve_x() FixNH::nve_x(); // update mu for dipoles - // d_mu/dt = omega cross mu - // renormalize mu to dipole length - + if (dipole_flag) { - double msq,scale,g[3]; double **mu = atom->mu; double **omega = atom->omega; int *mask = atom->mask; int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - if (mu[i][3] > 0.0) { - g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); - g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); - g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); - msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2]; - scale = mu[i][3]/sqrt(msq); - mu[i][0] = g[0]*scale; - mu[i][1] = g[1]*scale; - mu[i][2] = g[2]*scale; + if (dlm_flag == 0){ + // d_mu/dt = omega cross mu + // renormalize mu to dipole length + double msq,scale,g[3]; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (mu[i][3] > 0.0) { + g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); + g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); + g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); + msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2]; + scale = mu[i][3]/sqrt(msq); + mu[i][0] = g[0]*scale; + mu[i][1] = g[1]*scale; + mu[i][2] = g[2]*scale; + } + } else { + // Integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme + vector w, w_temp, a; + matrix Q, Q_temp, R; + double scale,s2,inv_len_mu; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit && mu[i][3] > 0.0) { + + // Construct Q from dipole: + // Q is the rotation matrix from space frame to body frame + // i.e. v_b = Q.v_s + + // Define mu to lie along the z axis in the body frame + // We take the unit dipole to avoid getting a scaling matrix + inv_len_mu = 1.0/mu[i][3]; + a[0] = mu[i][0]*inv_len_mu; + a[1] = mu[i][1]*inv_len_mu; + a[2] = mu[i][2]*inv_len_mu; + + // v = a x [0 0 1] - cross product of mu in space and body frames + // s = |v| + // c = a.[0 0 1] = a[2] + // vx = [ 0 -v[2] v[1] + // v[2] 0 -v[0] + // -v[1] v[0] 0 ] + // then + // Q = I + vx + vx^2 * (1-c)/s^2 + + s2 = a[0]*a[0] + a[1]*a[1]; + if (s2 != 0.0){ // i.e. the vectors are not parallel + scale = (1.0 - a[2])/s2; + + Q[0][0] = 1.0 - scale*a[0]*a[0]; Q[0][1] = -scale*a[0]*a[1]; Q[0][2] = -a[0]; + Q[1][0] = -scale*a[0]*a[1]; Q[1][1] = 1.0 - scale*a[1]*a[1]; Q[1][2] = -a[1]; + Q[2][0] = a[0]; Q[2][1] = a[1]; Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]); + } else { // if parallel then we just have I or -I + Q[0][0] = 1.0/a[2]; Q[0][1] = 0.0; Q[0][2] = 0.0; + Q[1][0] = 0.0; Q[1][1] = 1.0/a[2]; Q[1][2] = 0.0; + Q[2][0] = 0.0; Q[2][1] = 0.0; Q[2][2] = 1.0/a[2]; + } + + // Local copy of this particle's angular velocity (in space frame) + w[0] = omega[i][0]; w[1] = omega[i][1]; w[2] = omega[i][2]; + + // Transform omega into body frame: w_temp= Q.w + matvec(Q,w,w_temp); + + // Construct rotation R1 + BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]); + + // Apply R1 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R1 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Construct rotation R2 + BuildRyMatrix(R, dtf/force->ftm2v*w[1]); + + // Apply R2 to w: w_temp = R.w + matvec(R,w,w_temp); + + // Apply R2 to Q: Q = R^T.Q_temp + transpose_times3(R,Q_temp,Q); + + // Construct rotation R3 + BuildRzMatrix(R, 2.0*dtf/force->ftm2v*w_temp[2]); + + // Apply R3 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R3 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Construct rotation R4 + BuildRyMatrix(R, dtf/force->ftm2v*w[1]); + + // Apply R4 to w: w_temp = R.w + matvec(R,w,w_temp); + + // Apply R4 to Q: Q = R^T.Q_temp + transpose_times3(R,Q_temp,Q); + + // Construct rotation R5 + BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]); + + // Apply R5 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R5 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Transform w back into space frame w_temp = Q^T.w + transpose_matvec(Q_temp,w,w_temp); + omega[i][0] = w_temp[0]; omega[i][1] = w_temp[1]; omega[i][2] = w_temp[2]; + + // Set dipole according to updated Q: mu = Q^T.[0 0 1] * |mu| + mu[i][0] = Q_temp[2][0] * mu[i][3]; + mu[i][1] = Q_temp[2][1] * mu[i][3]; + mu[i][2] = Q_temp[2][2] * mu[i][3]; } + } + } } } diff --git a/src/fix_nve_sphere.cpp b/src/fix_nve_sphere.cpp index 3049605531..9f7b4a9eaa 100644 --- a/src/fix_nve_sphere.cpp +++ b/src/fix_nve_sphere.cpp @@ -21,13 +21,17 @@ #include "respa.h" #include "force.h" #include "error.h" +#include "math_vector.h" +#include "math_extra.h" using namespace LAMMPS_NS; using namespace FixConst; +using namespace MathExtra; #define INERTIA 0.4 // moment of inertia prefactor for sphere enum{NONE,DIPOLE}; +enum{NODLM,DLM}; /* ---------------------------------------------------------------------- */ @@ -41,13 +45,17 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) : // process extra keywords extra = NONE; + dlm = NODLM; int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"update") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nve/sphere command"); if (strcmp(arg[iarg+1],"dipole") == 0) extra = DIPOLE; - else error->all(FLERR,"Illegal fix nve/sphere command"); + else if (strcmp(arg[iarg+1],"dipole/dlm") == 0) { + extra = DIPOLE; + dlm = DLM; + } else error->all(FLERR,"Illegal fix nve/sphere command"); iarg += 2; } else error->all(FLERR,"Illegal fix nve/sphere command"); } @@ -57,7 +65,7 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) : if (!atom->sphere_flag) error->all(FLERR,"Fix nve/sphere requires atom style sphere"); if (extra == DIPOLE && !atom->mu_flag) - error->all(FLERR,"Fix nve/sphere dipole requires atom attribute mu"); + error->all(FLERR,"Fix nve/sphere update dipole requires atom attribute mu"); } /* ---------------------------------------------------------------------- */ @@ -83,8 +91,10 @@ void FixNVESphere::init() void FixNVESphere::initial_integrate(int vflag) { - double dtfm,dtirotate,msq,scale; + double dtfm,dtirotate,msq,scale,s2,inv_len_mu; double g[3]; + vector w, w_temp, a; + matrix Q, Q_temp, R; double **x = atom->x; double **v = atom->v; @@ -120,25 +130,127 @@ void FixNVESphere::initial_integrate(int vflag) omega[i][2] += dtirotate * torque[i][2]; } } - + // update mu for dipoles - // d_mu/dt = omega cross mu - // renormalize mu to dipole length + if (extra == DIPOLE) { double **mu = atom->mu; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - if (mu[i][3] > 0.0) { - g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); - g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); - g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); - msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2]; - scale = mu[i][3]/sqrt(msq); - mu[i][0] = g[0]*scale; - mu[i][1] = g[1]*scale; - mu[i][2] = g[2]*scale; + if (dlm == NODLM) { + // d_mu/dt = omega cross mu + // renormalize mu to dipole length + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (mu[i][3] > 0.0) { + g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]); + g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]); + g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]); + msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2]; + scale = mu[i][3]/sqrt(msq); + mu[i][0] = g[0]*scale; + mu[i][1] = g[1]*scale; + mu[i][2] = g[2]*scale; + } + } else { + // Integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit && mu[i][3] > 0.0) { + + // Construct Q from dipole: + // Q is the rotation matrix from space frame to body frame + // i.e. v_b = Q.v_s + + // Define mu to lie along the z axis in the body frame + // We take the unit dipole to avoid getting a scaling matrix + inv_len_mu = 1.0/mu[i][3]; + a[0] = mu[i][0]*inv_len_mu; + a[1] = mu[i][1]*inv_len_mu; + a[2] = mu[i][2]*inv_len_mu; + + // v = a x [0 0 1] - cross product of mu in space and body frames + // s = |v| + // c = a.[0 0 1] = a[2] + // vx = [ 0 -v[2] v[1] + // v[2] 0 -v[0] + // -v[1] v[0] 0 ] + // then + // Q = I + vx + vx^2 * (1-c)/s^2 + + s2 = a[0]*a[0] + a[1]*a[1]; + if (s2 != 0.0){ // i.e. the vectors are not parallel + scale = (1.0 - a[2])/s2; + + Q[0][0] = 1.0 - scale*a[0]*a[0]; Q[0][1] = -scale*a[0]*a[1]; Q[0][2] = -a[0]; + Q[1][0] = -scale*a[0]*a[1]; Q[1][1] = 1.0 - scale*a[1]*a[1]; Q[1][2] = -a[1]; + Q[2][0] = a[0]; Q[2][1] = a[1]; Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]); + } else { // if parallel then we just have I or -I + Q[0][0] = 1.0/a[2]; Q[0][1] = 0.0; Q[0][2] = 0.0; + Q[1][0] = 0.0; Q[1][1] = 1.0/a[2]; Q[1][2] = 0.0; + Q[2][0] = 0.0; Q[2][1] = 0.0; Q[2][2] = 1.0/a[2]; + } + + // Local copy of this particle's angular velocity (in space frame) + w[0] = omega[i][0]; w[1] = omega[i][1]; w[2] = omega[i][2]; + + // Transform omega into body frame: w_temp= Q.w + matvec(Q,w,w_temp); + + // Construct rotation R1 + BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]); + + // Apply R1 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R1 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Construct rotation R2 + BuildRyMatrix(R, dtf/force->ftm2v*w[1]); + + // Apply R2 to w: w_temp = R.w + matvec(R,w,w_temp); + + // Apply R2 to Q: Q = R^T.Q_temp + transpose_times3(R,Q_temp,Q); + + // Construct rotation R3 + BuildRzMatrix(R, 2.0*dtf/force->ftm2v*w_temp[2]); + + // Apply R3 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R3 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Construct rotation R4 + BuildRyMatrix(R, dtf/force->ftm2v*w[1]); + + // Apply R4 to w: w_temp = R.w + matvec(R,w,w_temp); + + // Apply R4 to Q: Q = R^T.Q_temp + transpose_times3(R,Q_temp,Q); + + // Construct rotation R5 + BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]); + + // Apply R5 to w: w = R.w_temp + matvec(R,w_temp,w); + + // Apply R5 to Q: Q_temp = R^T.Q + transpose_times3(R,Q,Q_temp); + + // Transform w back into space frame w_temp = Q^T.w + transpose_matvec(Q_temp,w,w_temp); + omega[i][0] = w_temp[0]; omega[i][1] = w_temp[1]; omega[i][2] = w_temp[2]; + + // Set dipole according to updated Q: mu = Q^T.[0 0 1] * |mu| + mu[i][0] = Q_temp[2][0] * mu[i][3]; + mu[i][1] = Q_temp[2][1] * mu[i][3]; + mu[i][2] = Q_temp[2][2] * mu[i][3]; } + } + } } } @@ -165,6 +277,7 @@ void FixNVESphere::final_integrate() // update v,omega for all particles // d_omega/dt = torque / inertia + double rke = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { dtfm = dtf / rmass[i]; @@ -176,5 +289,7 @@ void FixNVESphere::final_integrate() omega[i][0] += dtirotate * torque[i][0]; omega[i][1] += dtirotate * torque[i][1]; omega[i][2] += dtirotate * torque[i][2]; + rke += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2])*radius[i]*radius[i]*rmass[i]; } + } diff --git a/src/fix_nve_sphere.h b/src/fix_nve_sphere.h index 10a22be702..801518ff1d 100644 --- a/src/fix_nve_sphere.h +++ b/src/fix_nve_sphere.h @@ -34,6 +34,7 @@ class FixNVESphere : public FixNVE { protected: int extra; + int dlm; }; } @@ -53,12 +54,17 @@ E: Fix nve/sphere requires atom style sphere Self-explanatory. -E: Fix nve/sphere dipole requires atom attribute mu +E: Fix nve/sphere update dipole requires atom attribute mu An atom style with this attribute is needed. E: Fix nve/sphere requires extended particles This fix can only be used for particles of a finite size. + +E: Fix nve/sphere dlm must be used with update dipole + +The DLM algorithm can only be used in conjunction with update dipole. + */ diff --git a/src/force.cpp b/src/force.cpp index aa4ad817ad..37ac4f4ee6 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -260,12 +260,6 @@ char *Force::pair_match_ptr(Pair *ptr) if (ptr == hybrid->styles[i]) return hybrid->keywords[i]; } - if (strstr(pair_style,"hybrid/overlay")) { - PairHybridOverlay *hybrid = (PairHybridOverlay *) pair; - for (int i = 0; i < hybrid->nstyles; i++) - if (ptr == hybrid->styles[i]) return hybrid->keywords[i]; - } - return NULL; } diff --git a/src/math_extra.cpp b/src/math_extra.cpp index 7095412b3e..de61c8901c 100644 --- a/src/math_extra.cpp +++ b/src/math_extra.cpp @@ -607,6 +607,51 @@ void inertia_triangle(double *idiag, double *quat, double mass, inertia[5] = tensor[0][1]; } +/* ---------------------------------------------------------------------- + Build rotation matrix for a small angle rotation around the X axis + ------------------------------------------------------------------------- */ + +void BuildRxMatrix(double R[3][3], const double angle) +{ + const double angleSq = angle * angle; + const double cosAngle = (1.0 - angleSq * 0.25) / (1.0 + angleSq * 0.25); + const double sinAngle = angle / (1.0 + angleSq * 0.25); + + R[0][0] = 1.0; R[0][1] = 0.0; R[0][2] = 0.0; + R[1][0] = 0.0; R[1][1] = cosAngle; R[1][2] = -sinAngle; + R[2][0] = 0.0; R[2][1] = sinAngle; R[2][2] = cosAngle; +} + +/* ---------------------------------------------------------------------- + Build rotation matrix for a small angle rotation around the Y axis + ------------------------------------------------------------------------- */ + +void BuildRyMatrix(double R[3][3], const double angle) +{ + const double angleSq = angle * angle; + const double cosAngle = (1.0 - angleSq * 0.25) / (1.0 + angleSq * 0.25); + const double sinAngle = angle / (1.0 + angleSq * 0.25); + + R[0][0] = cosAngle; R[0][1] = 0.0; R[0][2] = sinAngle; + R[1][0] = 0.0; R[1][1] = 1.0; R[1][2] = 0.0; + R[2][0] = -sinAngle; R[2][1] = 0.0; R[2][2] = cosAngle; +} + +/* ---------------------------------------------------------------------- + Build rotation matrix for a small angle rotation around the Y axis + ------------------------------------------------------------------------- */ + +void BuildRzMatrix(double R[3][3], const double angle) +{ + const double angleSq = angle * angle; + const double cosAngle = (1.0 - angleSq * 0.25) / (1.0 + angleSq * 0.25); + const double sinAngle = angle / (1.0 + angleSq * 0.25); + + R[0][0] = cosAngle; R[0][1] = -sinAngle; R[0][2] = 0.0; + R[1][0] = sinAngle; R[1][1] = cosAngle; R[1][2] = 0.0; + R[2][0] = 0.0; R[2][1] = 0.0; R[2][2] = 1.0; +} + /* ---------------------------------------------------------------------- */ } diff --git a/src/math_extra.h b/src/math_extra.h index 7f0db9079f..4e6ee4c77c 100755 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -111,6 +111,10 @@ namespace MathExtra { inline void rotation_generator_x(const double m[3][3], double ans[3][3]); inline void rotation_generator_y(const double m[3][3], double ans[3][3]); inline void rotation_generator_z(const double m[3][3], double ans[3][3]); + + void BuildRxMatrix(double R[3][3], const double angle); + void BuildRyMatrix(double R[3][3], const double angle); + void BuildRzMatrix(double R[3][3], const double angle); // moment of inertia operations diff --git a/src/thermo.cpp b/src/thermo.cpp index 5127b91a89..5ecc161c8f 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -991,7 +991,7 @@ int Thermo::add_variable(const char *id) } /* ---------------------------------------------------------------------- - compute a single thermodyanmic value, word is any keyword in custom list + compute a single thermodynamic value, word is any keyword in custom list called when a variable is evaluated by Variable class return value as double in answer return 0 if str is recoginzed keyword, 1 if unrecognized