# Variable formula behaving erroneously

I am not able to see through the problem. My variable is defined as

``````variable wellFx atom -1*(sqrt(x*x+y*y)<=18)*(sqrt(x*x+y*y)>=10)*(2*\${k}*(sqrt(x*x+y*y)-10)*(-\${a0}+sqrt(x*x+y*y)-10)*(-\${a0}-\${a1}+sqrt(x*x+y*y)-10)*(x)*(1/(sqrt(x*x+y*y))))
``````

As it is defined here logically that `1/sqrt(x^2+y^2)` will be evaluated if and only `sqrt(x^2+y^2)` <18 or >10, still I get the error that variable attempted division by 0. How come it is possible. Are there any preferences in evaluating variable formula. I am working with phantom chain and there is an annular region where I have put a potential barrier. This formula works ok with self-avoiding chain but does not works with phantom chains.

Where does it say that?

LAMMPS does what the documentation says it does. So please provide proof that the behavior you see is in conflict with the documentation.

The logical expression `(sqrt(x*x+y*y)<=18)*(sqrt(x*x+y*y)>=10)` if `sqrt(x*x+y*y)` doesnâ€™t lie in this range??? Its evaluation will be 0 and the expression should then become 0. Isnâ€™t it? The behavior was fine for polymer with excluded-volume, but now it is giving error with a phantom polymer!

I am pasting below the script. I also tried using region command for the fix addforce, but still I get the same error. The region I define is an annular area with inner radius of 10 and outer radius of 18. Therefore, `sqrt(x*x+y*y)` can never be 0.

``````# VARIABLES
variable seed index 0
variable rand equal ceil(random(1423,8674437,\${seed}))
variable gamma_t string 3.0
variable gamma_r string 1.0
variable temp string 1.0
variable fp string 0.0
variable a0 equal 3
variable a02 equal 9
variable a1 equal 5.0
variable k equal 1.0
variable iens index 0
variable fens index 0
label start
variable ens loop \${iens} \${fens}
variable fname index init.txt

# Initialization
units		lj
dimension	2
boundary	p p p
atom_style	hybrid molecular sphere dipole
newton		off
log 		log.\${ens}.txt

# Potential information
neighbor	1.0 bin
neigh_modify	every 1 one 10000
bond_style      harmonic
bond_coeff	1 100 1.0

velocity 	all create \${temp} \${rand}

region          stable sphere 0.0 0.0 0.0 18 side out
region          metastable sphere 0.0 0.0 0.0 10.0 side in
region          inv_metastable  sphere 0.0 0.0 0.0 10.0 side out
region          barrier sphere 0.0 0.0 0.0 18 side in
region          shell_barrier intersect 2 inv_metastable barrier

#randomize dipole
set             type  * dipole/random \${rand} 1.0
write_dump	all custom writedump\${ens}.comp.cfg id x y z mux muy muz

group           atoms_equib id  40:260 1:31

timestep	0.0002

#force due to external pot
variable wellFx atom -1*(sqrt(x*x+y*y)<=18)*(sqrt(x*x+y*y)>=10)*(2*\${k}*(sqrt(x*x+y*y)-10)*(-\${a0}+sqrt(x*x+y*y)-10)*(-\${a0}-\${a1}+sqrt(x*x+y*y)-10)*(x)*(1/(sqrt(x*x+y*y)))) #-0.01*(z-x))
variable wellFy atom -1*(sqrt(x*x+y*y)<=18)*(sqrt(x*x+y*y)>=10)*(2*\${k}*(sqrt(x*x+y*y)-10)*(-\${a0}+sqrt(x*x+y*y)-10)*(-\${a0}-\${a1}+sqrt(x*x+y*y)-10)*(y)*(1/(sqrt(x*x+y*y)))) #-0.01*(z-x))
variable boundy atom -1*(x>=10.0)*(x<=18)*(y<=0.5)*(y>=-0.5)*(10000.0*(y)) #-0.01*(z-x))

fix         step atoms_equib brownian/sphere \${temp} \${rand} rng gaussian gamma_t \${gamma_t} gamma_r \${gamma_r}
# self-propulsion force along the dipole direction
#fix         activity all propel/self dipole \${fp}
fix             momentum all momentum 100 linear 1 1 0
fix             2d all enforce2d
thermo_style    custom atoms step temp
thermo          1000
dump            equib1 all atom 5000 equib\${ens}.lammpstrj
dump_modify     equib1 time yes
run             100000
unfix step
#unfix activity
unfix momentum
unfix 2d
undump equib1
``````

Please point to the section in the documentation for the variable command that its evaluation of variable expression behaves like this.

What you describe is behavior that is not necessarily true for any programming language. For example it is true for C/C++ but not for Fortran.

Well the error occurs in variable.cpp then it should work! The behavior is true to C++ for excluded volume chains. But suddenly for phantom chain it gives this error.

You are missing the point here. While the LAMMPS variable syntax is implemented in C++ it does not automatically mean that it follows the C++ short circuiting rules.

In this script, I am using region argument for fix addforce. That should avoid the question of erroneous behavior for the variable definition. But still I get the same error. Division by zero, but no particles in this region can have radius equal to zero. Then how come it gives the division by zero error.

``````variable seed index 0
variable rand equal ceil(random(1423,8674437,\${seed}))
variable gamma_t string 3.0
variable gamma_r string 1.0
variable temp string 1.0
variable fp string 0.0
variable a0 equal 3
variable a02 equal 9
variable a1 equal 5.0
variable k equal 1.0
variable iens index 0
variable fens index 0
label start
variable ens loop \${iens} \${fens}
variable fname index init.txt

# Initialization
units		lj
dimension	2
boundary	p p p
atom_style	hybrid molecular sphere dipole
newton		off
log 		log.\${ens}.txt

# Potential information
neighbor	1.0 bin
neigh_modify	every 1 one 10000
bond_style      harmonic
bond_coeff	1 100 1.0

velocity 	all create \${temp} \${rand}

region          stable sphere 0.0 0.0 0.0 18 side out
region          metastable sphere 0.0 0.0 0.0 10.0 side in
region          inv_metastable  sphere 0.0 0.0 0.0 10.0 side out
region          barrier sphere 0.0 0.0 0.0 18 side in
region          shell_barrier intersect 2 inv_metastable barrier

#randomize dipole
set             type  * dipole/random \${rand} 1.0
write_dump	all custom writedump\${ens}.comp.cfg id x y z mux muy muz

group           atoms_equib id  40:260 1:31

timestep	0.0002

#force due to external pot
variable wellFx atom -1*(2*\${k}*(sqrt(x*x+y*y)-10)*(-\${a0}+sqrt(x*x+y*y)-10)*(-\${a0}-\${a1}+sqrt(x*x+y*y)-10)*x*(1/(sqrt(x*x+y*y)))) #-0.01*(z-x))
variable wellFy atom -1*(2*\${k}*(sqrt(x*x+y*y)-10)*(-\${a0}+sqrt(x*x+y*y)-10)*(-\${a0}-\${a1}+sqrt(x*x+y*y)-10)*y*(1/(sqrt(x*x+y*y)))) #-0.01*(z-x))
variable boundy atom -1*(x>=10.0)*(x<=18)*(y<=0.5)*(y>=-0.5)*(10000.0*(y)) #-0.01*(z-x))

fix         step atoms_equib brownian/sphere \${temp} \${rand} rng gaussian gamma_t \${gamma_t} gamma_r \${gamma_r}
# self-propulsion force along the dipole direction
#fix         activity all propel/self dipole \${fp}
fix             momentum all momentum 100 linear 1 1 0
fix             2d all enforce2d
thermo_style    custom atoms step temp
thermo          1000
dump            equib1 all atom 5000 equib\${ens}.lammpstrj
dump_modify     equib1 time yes
run             100000
unfix step
#unfix activity
unfix momentum
unfix 2d
undump equib1

``````

I get your point, thatâ€™s why I am adding now the region argument for fix addforce, but still get the same error. The region clearly defines with minimum radius of 10. No way a particle in the region can have 0 radius.

init.txt (10.7 KB)
The input topology for the polymer.

Again, this comes down to the documentation. Does the documentation for fix addforce explicitly state anywhere that it will evaluate variables after the check whether a particle is within the region and only for those atoms?

From an implementation point of view it makes rather little sense, since atom style variable evaluations are only efficient if done for all atoms.

You are arguing with common sense, but that does not apply to computer programs, they strictly follow what is in the program.

My common sense is working all good for self-avoiding polymer. But, with the same setup if I just remove pair style and put harmonic bonds instead of FENE (with same equilibrium bond length), I get the error. If there was any difference in common sense and the algorithm, then it should fail every time, and one would be always in doubt whether the algorithm will follow the logic in any simulation setup! I am pretty sure the algorithms are written to follow the logic.

Yes, the documentation clearly says that the force will be added only to atoms present in the region.

If the region keyword is used, the atom must also be in the specified geometric region in order to have force added to it.

Then obviously there must be a check before adding force to an atom, whether it belongs to the region or not!

The documentation says that the force is applied only to atoms in the region. But there is no mention of the force being only computed for the atoms in the region.

In fact, there is no function or feature in the variable code to have atom style variables only evaluated for atoms in specific groups or regions.

You are arguing common sense and not computer logic again.

Let us leave aside this region thing. I understand there are several executions involved.
But, I want to understand how can this expression go wrong?

``````variable wellFx atom -1*(sqrt(x*x+y*y)<=18)*(sqrt(x*x+y*y)>=10)*(2*\${k}*(sqrt(x*x+y*y)-10)*(-\${a0}+sqrt(x*x+y*y)-10)*(-\${a0}-\${a1}+sqrt(x*x+y*y)-10)*(x)*(1/(sqrt(x*x+y*y))))
``````

I used atom style variable, and the logical expressions of ` (sqrt(x*x+y*y)<=18) and (sqrt(x*x+y*y)>=10)` is computed for each atomâ€™s positions, and will return 0 or 1. Then should I do not expect this composite expression to provide me a correct solution? In case the radius is 0, the resultant solution to the expression should be 0 and not the error?

If you want somebody to have a closer look you have to do the following:

• name the LAMMPS version you are using, the platform you are running on, how compiled, which command line flags etc. the stuff that should have been at the beginning of your very first post. Please re-read the forum guidelines (completely!). We put them there for a reason. If you want to be treated nicely and want to get competent help, then you have to play by those guidelines. Otherwise, you are likely to be ignored.
• create a minimum input deck that reports only the error you are concerned about. When I run the files you posted, I get:
``````./lmp -in in.atomstyle
LAMMPS (27 Jun 2024)
WARNING: Atom style hybrid defines both, per-type and per-atom masses; both must be set, but only per-atom masses will be used (src/atom_vec_hybrid.cpp:132)
orthogonal box = (-100 -50 -1) to (100 50 1)
1 by 1 by 1 MPI processor grid
260 atoms
scanning bonds ...
2 = max bonds/atom
orthogonal box = (-100 -50 -1) to (100 50 1)
1 by 1 by 1 MPI processor grid
259 bonds
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj:    0        0        0
special bond factors coul:  0        0        0
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000 seconds
ERROR: Variable rand: Invalid math function in variable formula (src/variable.cpp:3758)
Last command: velocity 	all create \${temp} \${rand}
``````

You are reporting a division by zero error and in your variable expression this can happen when both x and y are zero for an atom since you have:

Since you are creating atoms on a regular unshifted lattice, this is quite likely to happen,
so I donâ€™t see that there is much of a mystery at all, but just making assumptions that are not backed up by the documentation (like assuming short-circuiting on variable expressions).

I have now repeatedly asked for you to back up your claims by pointing to confirmation in the documentation and you have failed to do so.

For some molecules, sqrt(x^2+y^2) may be equal to 0. My approach is to use 1/sqrt(x^2+y^2+1E-20) to avoid this low probability event. I hope this solution is helpful for your issue.

1 Like

Yes, I have already added a new variable which takes care of zero for other atoms not in the region.