ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

change to 9 element stress tensor

Browse Source
This commit is contained in:
Steve Plimpton
2022-06-16 17:38:57 -06:00
parent 4dbecbba51
commit 1f35065afc
4 changed files with 24 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
examples/mdi/sequence_driver.py
View File

@ -155,7 +155,7 @@ def perform_tasks(world,mdicomm,dummy):
# request virial tensor
mdi.MDI_Send_command("<STRESS",mdicomm)
virial = mdi.MDI_Recv(6,mdi.MDI_DOUBLE,mdicomm)
virial = mdi.MDI_Recv(9,mdi.MDI_DOUBLE,mdicomm)
virial = world.bcast(virial,root=0)
# request forces
@ -165,10 +165,10 @@ def perform_tasks(world,mdicomm,dummy):
world.Bcast(forces,root=0)
# final output from each calculation
# pressure = just virial component, no kinetic component
# pressure = trace of virial tensor, no kinetic component
aveeng = pe/natoms
pressure = (virial[0] + virial[1] + virial[2]) / 3.0
pressure = (virial[0] + virial[4] + virial[8]) / 3.0
m = 0
fx = fy = fz = 0.0

18
src/MDI/fix_mdi_qm.cpp
View File

@ -329,18 +329,22 @@ void FixMDIQM::post_force(int vflag)
}
}
// optionally request pressure tensor from MDI engine, convert to virial
// optionally request stress tensor from MDI engine, convert to virial
// qm_virial = fix output for global QM virial
if (virialflag) {
ierr = MDI_Send_command("<STRESS", mdicomm);
if (ierr) error->all(FLERR, "MDI: <STRESS command");
ierr = MDI_Recv(qm_virial, 6, MDI_DOUBLE, mdicomm);
ierr = MDI_Recv(qm_virial, 9, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: <STRESS data");
MPI_Bcast(qm_virial, 6, MPI_DOUBLE, 0, world);
MPI_Bcast(qm_virial, 9, MPI_DOUBLE, 0, world);
for (int i = 0; i < 6; i++)
qm_virial[i] *= mdi2lmp_pressure;
qm_virial_symmetric[0] = qm_virial[0] * mdi2lmp_pressure;
qm_virial_symmetric[1] = qm_virial[4] * mdi2lmp_pressure;
qm_virial_symmetric[2] = qm_virial[8] * mdi2lmp_pressure;
qm_virial_symmetric[3] = 0.5*(qm_virial[1]+qm_virial[3]) * mdi2lmp_pressure;
qm_virial_symmetric[4] = 0.5*(qm_virial[2]+qm_virial[6]) * mdi2lmp_pressure;
qm_virial_symmetric[5] = 0.5*(qm_virial[5]+qm_virial[7]) * mdi2lmp_pressure;
}
// optionally set fix->virial
@ -356,7 +360,7 @@ void FixMDIQM::post_force(int vflag)
else if (domain->dimension == 3)
volume = domain->xprd * domain->yprd * domain->zprd;
for (int i = 0; i < 6; i++)
virial[i] = qm_virial[i]*volume/nprocs;
virial[i] = qm_virial_symmetric[i]*volume/nprocs;
}
}
@ -382,7 +386,7 @@ double FixMDIQM::compute_scalar()
double FixMDIQM::compute_vector(int n)
{
return qm_virial[n];
return qm_virial_symmetric[n];
}
/* ----------------------------------------------------------------------

2
src/MDI/fix_mdi_qm.h
View File

@ -47,7 +47,7 @@ class FixMDIQM : public Fix {
double qm_energy;
int lmpunits;
double qm_virial[6],qm_virial_all[6];
double qm_virial[9],qm_virial_symmetric[6];
double **fqm;
MDI_Comm mdicomm;

13
src/MDI/mdi_engine.cpp
View File

@ -1413,16 +1413,21 @@ void MDIEngine::send_pe()
/* ----------------------------------------------------------------------
<STRESS command
send 6-component stress tensor (no kinetic energy term)
send 9-component stress tensor (no kinetic energy term)
---------------------------------------------------------------------- */
void MDIEngine::send_stress()
{
double vtensor[6];
double vtensor_full[9];
press->compute_vector();
for (int i = 0; i < 6; i++) vtensor[i] = press->vector[i] * lmp2mdi_pressure;
vtensor_full[0] = press->vector[0] * lmp2mdi_pressure;
vtensor_full[4] = press->vector[1] * lmp2mdi_pressure;
vtensor_full[8] = press->vector[2] * lmp2mdi_pressure;
vtensor_full[1] = vtensor_full[3] = press->vector[3] * lmp2mdi_pressure;
vtensor_full[2] = vtensor_full[6] = press->vector[4] * lmp2mdi_pressure;
vtensor_full[5] = vtensor_full[7] = press->vector[5] * lmp2mdi_pressure;
int ierr = MDI_Send(vtensor, 6, MDI_DOUBLE, mdicomm);
int ierr = MDI_Send(vtensor_full, 9, MDI_DOUBLE, mdicomm);
if (ierr) error->all(FLERR, "MDI: <STRESS data");
}