Scaling of Langevin Drag force when using metal units

LAMMPS LAMMPS Beginners
1

Dear All,

The drag force applied by fix langevin at T=0 when using metal units is different compared to LJ units. I am unsure how this fits with the definition from the ‘fix langevin’ page where the drag force is F_d=-m*v/g, details below:

Setup:

  • 3D domain with 1 atom initialized with velocity (vx,vy,vz)=(0,0,-1) (ps/angstrom or dimensionless in LJ case)
  • I apply: fix Langevin with T=0, damping factor, g=1 (ps or dimensionless in LJ case).
  • Use fix nve and look at the force at time t=0.

Output:

  1. Case 1: When using LJ units with m=1, g=1, I get a force of (fx,fy,fx)=(0,0,1). This matches the expected drag force of F_d=-m*v/g, from the ‘fix langevin’ documentation
  2. Case 2: When using metal units, m=1 (g/mole) , g=1(ps) I get a force of (fx,fy,fx)=(0,0,0.000103643)(ev/Angstrom), I am not sure where this additional factor of 1.03x10^-4 comes from? This may be just my misunderstanding of how to convert these units e.g. how many moles are contained within each sphere which have diameter 1angstrom?

version & platform: LAMMPS (2 Aug 2023), solved using lmp_serial on Apple M2 processor.

Thank you for any help!

Case 2 File (Metal Units)

units metal
atom_style sphere

lattice fcc 3
region mybox block 0 5 0 5 0 5 units lattice
create_box 1 mybox
create_atoms 1 single 1 1 1

set atom * mass 1
set atom * vx 0.0 vy 0.0 vz -1.0

fix 1 all langevin 0 0 1 8927 omega no # apply langevin damping
fix 2 all nve

timestep 0.001

compute 1 all property/atom fx fy fz radius mass 

dump 1 all custom 1 dump.solution id x y z vx vy vz c_1[*]

run 8
2

You have to look at the source code to work this out.

You have:

        fdrag[0] = gamma1*v[i][0];
        fdrag[1] = gamma1*v[i][1];
        fdrag[2] = gamma1*v[i][2];

With:

     gamma1 = gfactor1[type[i]];

And:

   gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;

The ftm2v variable is the unit specific scaling factor that converts a force*time/mass expression to velocity units and it is defined in the set_units() method in the Update class. For metal units you have:

    force->ftm2v = 1.0 / 1.0364269e-4;
3

Dear Axel,

Thank you for that fast reply!

Am I correct in thinking that this factor ftm2v will also appear with the random forces if a non zero temperature is used?

I had a look at the source code and it seems the random forces ‘fran’ are multiplied by ‘gfactor2[i]’
which is defined (Line 292) of fix_langevin.cpp

      gfactor2[i] = sqrt(atom->mass[i]) / force->ftm2v;

which is then used on lines 661-663:

 fran[0] = gamma2*(random->uniform()-0.5);
 fran[1] = gamma2*(random->uniform()-0.5);
 fran[2] = gamma2*(random->uniform()-0.5);

Thanks again!

Edwina