Dear all, especially dear Axel and Steve,

I have posted here previously about adding external one body potential V(x,y,z) for a system.

There was two suggestion and in the end I took Axel’s suggestion to implement it through

fix addforce and atom-style variables.

When I tried to test it I found that:

a) either I do not understand how to print out the potential energy of the system

b) or fix addforce does not add the energy to the potential energy for some reason

c) I have messed up something in the syntax of adding the potential to the system through “energy”

Please help me find out which is the case! Please note that when I print out the potential

it gives reasonable values, but potential energy with the “pe” variable always gives 0.

I do not ask you to check the details of the potential it is only provided for illustration.

Please find the relevant input bellow.

Thank you for your help in anticipation!

Janos

#input:

include “system.in.init”

#(…)

# PARAMETERS

variable pi equal 3.14159

variable b equal 4.83045

variable Ea equal 0.00698237

variable Ec equal 3.6749e-05

variable bb equal 0.423342

variable y0 equal 37.9585

#calculating the potential: PROBABLY NOT IMPORTANT ONLY FUNCTIONS OF x, y and z

variable Axz atom ({Ea}-{Ec}*(6.0-4.0*cos(2.0*${pi}*x/b)*cos(2.0*{pi}**z/b/sqrt(3.0))-2.0*cos(4.0*{pi}**z/b/sqrt(3.0)))/9.0)
variable By atom (exp(-2.0*{bb}*(y-{y0})-2.0*exp(-{bb}(y-${y0}))))

variable Vxyz atom (v_Axzv_By)

#calculating derivatives: PROBABLY NOT IMPORTANT ONLY FUNCTIONS OF x, y and z

variable CuFx atom (-1.0*{Ec}*8.0/9.0*{pi}/b*cos(2.0*{pi}*z/b/sqrt(3.0))*sin(2.0*{pi}**x/b)*v_By)
variable CuFy atom (-1.0*v_Axz*(-2.0*{bb}**exp(-2.0*{bb}*(y-{y0}))+2{bb}*exp(-1*{bb}(y-{y0}))))
variable CuFz atom (-1.0*{Ec}*8.0/9.0*${pi}/sqrt(3.0)/b*sin(4.0*{pi}*z/sqrt(3.0)/$b)*v_By)

#Apply the force on the Xe atom: PLEASE CHECK THIS PART

fix field Xe addforce v_CuFx v_CuFy v_CuFz energy v_Vxyz

#Printing out quantities: PLEASE CHECK THIS PART

compute ener all reduce sum v_Vxyz

variable e equal c_ener

variable p equal pe

thermo 1

thermo_style custom pe c_ener #this gives 0 and some meaningful number for c_ener

fix extra all print 1 “{p} {e}” file ener.dat #this one does the same

#one note:

#the system is set up with no other potential then the one body V(x,y,z)

#eg in system.in.init

units electron

atom_style full

neigh_modify exclude type 1 1

kspace_style none

atom_modify sort 1 1.0

pair_style none