diff --git a/examples/cij_nvt/in.ljfcc b/examples/cij_nvt/in.ljfcc index 9583f113a7..62c81d060f 100644 --- a/examples/cij_nvt/in.ljfcc +++ b/examples/cij_nvt/in.ljfcc @@ -11,7 +11,7 @@ variable nsteps index 50 # length of run variable nthermo index 10 # thermo output interval -variable nlat index 3 # size of box +variable nlat index 5 # size of box variable delta index 1.0e-6 # strain size variable temp index 0.3 # temperature variable rho index 0.8442 # reduced density @@ -38,10 +38,9 @@ fix 1 all nvt temp ${temp} ${temp} 0.1 run 100 compute born all born/matrix -variable born vector c_born/vol thermo ${nthermo} -thermo_style custom step temp pe press v_born[1] v_born[2] v_born[3] v_born[4] v_born[5] v_born[6] v_born[7] v_born[8] v_born[9] v_born[10] v_born[11] v_born[12] v_born[13] v_born[14] v_born[15] v_born[16] v_born[17] v_born[18] v_born[19] v_born[20] v_born[21] -thermo_modify line multi +thermo_style custom step temp pe press c_born[*] +thermo_modify line multi run ${nsteps} diff --git a/examples/cij_nvt/in.ljfcc_both b/examples/cij_nvt/in.ljfcc_both index 2cfc3107f5..f83484341f 100644 --- a/examples/cij_nvt/in.ljfcc_both +++ b/examples/cij_nvt/in.ljfcc_both @@ -13,7 +13,7 @@ variable nsteps index 200 # length of run variable nthermo index 10 # thermo output interval -variable nlat index 3 # size of box +variable nlat index 5 # size of box variable delta index 1.0e-8 # strain size variable temp index 0.3 # temperature variable rho index 0.8442 # reduced density @@ -37,15 +37,14 @@ neigh_modify every 1 delay 0 check no fix 1 all nvt temp ${temp} ${temp} 0.1 -variable vol equal vol compute myvirial all pressure NULL virial compute bornnum all born/matrix numdiff ${delta} myvirial compute born all born/matrix -variable bornrel vector 1.0-c_bornnum/c_born*v_vol -fix bornrel all ave/time 1 1 ${nthermo} v_bornrel mode vector ave running +variable bornrel vector 1.0-c_bornnum/c_born +fix bornrel all ave/time 1 1 ${nthermo} v_bornrel start 100 mode vector ave running thermo ${nthermo} -thermo_style custom step temp pe press & +thermo_style custom step temp pe press c_born[10] c_bornnum[10] & v_bornrel[1] v_bornrel[2] v_bornrel[3] v_bornrel[4] v_bornrel[5] v_bornrel[6] v_bornrel[7] v_bornrel[8] v_bornrel[9] v_bornrel[10] v_bornrel[11] v_bornrel[12] v_bornrel[13] v_bornrel[14] v_bornrel[15] v_bornrel[16] v_bornrel[17] v_bornrel[18] v_bornrel[19] v_bornrel[20] v_bornrel[21] thermo_modify line multi diff --git a/examples/cij_nvt/in.ljfcc_num b/examples/cij_nvt/in.ljfcc_num index c45f3aba80..a9ad8c41ae 100644 --- a/examples/cij_nvt/in.ljfcc_num +++ b/examples/cij_nvt/in.ljfcc_num @@ -13,7 +13,7 @@ variable nsteps index 50 # length of run variable nthermo index 10 # thermo output interval -variable nlat index 3 # size of box +variable nlat index 5 # size of box variable delta index 1.0e-6 # strain size variable temp index 0.3 # temperature variable rho index 0.8442 # reduced density diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp index b92eb5bd40..a1b30fc421 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.cpp +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -633,8 +633,8 @@ void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) // A.T. // this works for vector indices 7, 8, 9, 12, 14, 18 and 15, 16, 17 // corresponding i,j indices 12, 13, 14, 23, 25, 36 and 26, 34, 35 - int k = dirlist[idir][1]; - int l = dirlist[idir][0]; + // int k = dirlist[idir][1]; + // int l = dirlist[idir][0]; // A.T. // this works for vector indices 7, 8, 9, 12, 14, 18 and 10, 11, 13 @@ -642,8 +642,8 @@ void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) // G.C.: // I see no difference with a 0 step simulation between both // methods. - // int k = dirlist[idir][0]; - // int l = dirlist[idir][1]; + int k = dirlist[idir][0]; + int l = dirlist[idir][1]; for (int i = 0; i < nall; i++) x[i][k] = temp_x[i][k] + numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); } @@ -715,29 +715,61 @@ void ComputeBornMatrix::virial_addon() double* sigv = compute_virial->vector; double modefactor[6] = {1.0, 1.0, 1.0, 0.5, 0.5, 0.5}; - for (int idir = 0; idir < NDIR_VIRIAL; idir++) { - int ijvgt = idir; // this is it. - double addon; + // Back to the ugly way + // You can compute these factor by looking at + // every Dijkl terms and adding the proper virials + // Take into account the symmetries. For example: + // B2323 = s33+D2323; B3232= s22+D3232; + // but D3232=D2323 (computed in compute_numdiff) + // and Cijkl = (Bijkl+Bjikl+Bijlk+Bjilk)/4. = (Bijkl+Bjilk)/2. + // see Yoshimoto eq 15.and eq A3. + values_global[0] += 2.0*sigv[0]; + values_global[1] += 2.0*sigv[1]; + values_global[2] += 2.0*sigv[2]; + values_global[3] += 0.5*(sigv[1]+sigv[2]); + values_global[4] += 0.5*(sigv[0]+sigv[2]); + values_global[5] += 0.5*(sigv[0]+sigv[1]); + values_global[6] += 0.0; + values_global[7] += 0.0; + values_global[8] += 0.0; + values_global[9] += sigv[4]; + values_global[10] += sigv[3]; + values_global[11] += 0.0; + values_global[12] += sigv[5]; + values_global[13] += 0.0; + values_global[14] += sigv[3]; + values_global[15] += sigv[5]; + values_global[16] += sigv[4]; + values_global[17] += 0.0; + values_global[18] += 0.0; + values_global[19] += sigv[4]; + values_global[20] += sigv[5]; - // extract the two indices composing the voigt representation + // This loop is actually bogus. + // + // for (int idir = 0; idir < NDIR_VIRIAL; idir++) { + // int ijvgt = idir; // this is it. + // double addon; - id = voigt3VtoM[ijvgt][0]; - jd = voigt3VtoM[ijvgt][1]; + // // extract the two indices composing the voigt representation - for (int knvgt=ijvgt; knvgt < NDIR_VIRIAL; knvgt++) { - kd = voigt3VtoM[knvgt][0]; - nd = voigt3VtoM[knvgt][1]; - addon = kronecker[id][nd]*sigv[virialMtoV[jd][kd]] + - kronecker[id][kd]*sigv[virialMtoV[jd][nd]]; - if(id != jd) - addon += kronecker[jd][nd]*sigv[virialMtoV[id][kd]] + - kronecker[jd][kd]*sigv[virialMtoV[id][nd]]; + // id = voigt3VtoM[ijvgt][0]; + // jd = voigt3VtoM[ijvgt][1]; - m = revalbe[ijvgt][knvgt]; + // for (int knvgt=ijvgt; knvgt < NDIR_VIRIAL; knvgt++) { + // kd = voigt3VtoM[knvgt][0]; + // nd = voigt3VtoM[knvgt][1]; + // addon = kronecker[id][nd]*sigv[virialMtoV[jd][kd]] + + // kronecker[id][kd]*sigv[virialMtoV[jd][nd]]; + // if(id != jd) + // addon += kronecker[jd][nd]*sigv[virialMtoV[id][kd]] + + // kronecker[jd][kd]*sigv[virialMtoV[id][nd]]; - values_global[revalbe[ijvgt][knvgt]] += 0.5*modefactor[idir]*addon; - } - } + // m = revalbe[ijvgt][knvgt]; + + // values_global[revalbe[ijvgt][knvgt]] += 0.5*modefactor[idir]*addon; + // } + // } } /* ----------------------------------------------------------------------