The pressure tensor calculated from virial is different from that of ThermodynamicQuantities

[Work Ren]

Hi,

I am trying to calculate the pressure tensors from virials and compare that from ThermodynamicQuantities, and found they are always different, although there are correlation between them, but the absolute values are different. Is this related to GPU calculation? or anything I missed?

Results fromThermodynamicQuantities/pressure: 

 0.36219389239583655 0.36198863432314066 0.3616700666714205 0.3619508644634659 

Results from kinetics and virial: 0.346732142010379 0.3425870620841554 0.36186231589715995 0.35039383999723145

Thank you very much.

Best Regards

Ren Hua

I posted the whole script for your reference.

#____start____

import gsd.hoomd

import hoomd

import matplotlib

import numpy

import hoomd

import math

import os

import itertools

%matplotlib inline

matplotlib.style.use('ggplot')

hoomd.version.version

import datetime

class Status():

    def __init__(self, sim):

        self.sim = sim

    @property

    def seconds_remaining(self):

        try:

            return (self.sim.final_timestep - self.sim.timestep) / self.sim.tps

        except ZeroDivisionError:

            return 0

    @property

    def etr(self):

        return str(datetime.timedelta(seconds=self.seconds_remaining))

gpu = hoomd.device.GPU()

sim = hoomd.Simulation(device=gpu)

sim.create_state_from_gsd(filename='random.gsd')

integrator = hoomd.md.Integrator(dt=0.005)

cell = hoomd.md.nlist.Cell()

lj = hoomd.md.pair.LJ(nlist=cell)

lj.params[('A', 'A')] = dict(epsilon=1, sigma=1)

lj.r_cut[('A', 'A')] = 2.5

integrator.forces.append(lj)

nvt = hoomd.md.methods.NVT(kT=1.5, filter=hoomd.filter.All(), tau=1.0)

integrator.methods.append(nvt)

sim.operations.integrator = integrator

logger_table = hoomd.logging.Logger(categories=['scalar', 'string'])

logger_table.add(sim, quantities=['timestep', 'tps'])

status = Status(sim)

logger_table[('Status', 'etr')] = (status, 'etr', 'string')

table = hoomd.write.Table(trigger=hoomd.trigger.Periodic(period=10000),

                          logger=logger_table)

logger1 = hoomd.logging.Logger()

logger1.add(lj, quantities=['energies', 'forces','torques','virials'])

therm_properties = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())

sim.operations.computes.append(therm_properties)

logger1.add(therm_properties)

gsd_writer1 = hoomd.write.GSD(filename='trajectory.gsd',

                             trigger=hoomd.trigger.Periodic(1000),

                             mode='ab',

                             filter=hoomd.filter.All(),

                             log=logger1)

sim.operations.writers.append(gsd_writer1)

sim.operations.writers.append(table)

sim.run(1000000)

import numpy as np

traj = gsd.hoomd.open('trajectory.gsd', 'rb')

lx1=traj[-1].configuration.box[0]

ly1=traj[-1].configuration.box[1]

lz1=traj[-1].configuration.box[2]

volume_box = lx1 * ly1 * lz1

pressure_xx_therm =[]

pressure_yy_therm =[]

pressure_zz_therm =[]

pressure_therm =[]

for traj1 in traj:

    pressure_tensor1 = traj1.log['md/compute/ThermodynamicQuantities/pressure_tensor']

    pressure1 = traj1.log['md/compute/ThermodynamicQuantities/pressure']

    pressure_xx_therm.append(pressure_tensor1[0])

    pressure_yy_therm.append(pressure_tensor1[3])

    pressure_zz_therm.append(pressure_tensor1[5])

    

    pressure_therm.append(pressure1[0])

pressure_xx_therm = np.array(pressure_xx_therm)

pressure_yy_therm = np.array(pressure_yy_therm)

pressure_zz_therm = np.array(pressure_zz_therm)

pressure_therm = np.array(pressure_therm)

# calculate the pressure tensor from DPD virial

pressure_xx_virial =[]

pressure_yy_virial  =[]

pressure_zz_virial  =[]

pressure_virial  =[]

for traj1 in traj:

    virials_lj = traj1.log['particles/md/pair/LJ/virials']

    

    v_xx_lj = virials_lj[:, 0]

    v_yy_lj = virials_lj[:, 3]

    v_zz_lj = virials_lj[:, 5]

    velocity = traj1.particles.velocity

    p_mass = traj1.particles.mass

    velocity_x = traj1.particles.velocity[:,0]

    velocity_x_power2 = np.power(velocity_x,2)

    dynamic_x = np.multiply(p_mass, velocity_x_power2)

    velocity_y =traj1.particles.velocity[:,1]

    velocity_y_power2 = np.power(velocity_y,2)

    dynamic_y = np.multiply(p_mass, velocity_y_power2)

    velocity_z = traj1.particles.velocity[:,2]

    velocity_z_power2 = np.power(velocity_z,2)

    dynamic_z = np.multiply(p_mass, velocity_z_power2)

    total_x= np.sum(dynamic_x)  + np.sum(v_xx_lj) 

    total_y= np.sum(dynamic_y)  + np.sum(v_yy_lj) 

    total_z= np.sum(dynamic_z)  + np.sum(v_zz_lj) 

    p_x_cal = total_x/volume_box

    p_y_cal = total_y/volume_box

    p_z_cal = total_z/volume_box

    pressure_xx_virial.append(p_x_cal)

    pressure_yy_virial.append(p_y_cal)

    pressure_zz_virial.append(p_z_cal)

    pressure_virial.append(np.mean([p_x_cal, p_y_cal, p_z_cal]))

pressure_xx_virial = np.array(pressure_xx_virial)

pressure_yy_virial = np.array(pressure_yy_virial)

pressure_zz_virial = np.array(pressure_zz_virial)

pressure_virial = np.array(pressure_virial)

print('Results fromThermodynamicQuantities/pressure:\n', np.mean(pressure_xx_therm), np.mean(pressure_yy_therm), np.mean(pressure_zz_therm), np.mean(pressure_therm))

print('Results from kinetics and virial:\n', np.mean(pressure_xx_virial), np.mean(pressure_yy_virial), np.mean(pressure_zz_virial), np.mean(pressure_virial))

#____end_____

[Work Ren]

The image of comparison and correlation

[Joshua Anderson]

You need to pass the argument `dynamic=['momentum']` to hoomd.write.GSD.

By default, GSD saves space and does not write particle velocities out on every frame. Your loop over trajectory frames is re-using the velocities from the first frame.
------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> On Jun 14, 2021, at 8:00 AM, Work Ren <rh....@gmail.com> wrote:
>
> The image of comparison and correlation<2021-06-14 11-16-14屏幕截图.png>

> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/c3f2d30f-371b-415d-afc8-4e6fdcbb7587n%40googlegroups.com.
> <2021-06-14 11-16-14屏幕截图.png>

[Work Ren]

Hi Joshua,

Thank you very much for the quick and helpful response, this is exactly what I am looking for, very helpful, thank you.