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.
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
read_data ${fname}
# 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 addF1 all addforce v_wellFx v_wellFy 0.0 region inv_metastable
fix addF2 all addforce 0.0 v_boundy 0.0
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
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
read_data ${fname}
# 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 addF1 all addforce v_wellFx v_wellFy 0.0 region shell_barrier
fix addF2 all addforce 0.0 v_boundy 0.0
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.
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.
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)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
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)
Reading data file ...
orthogonal box = (-100 -50 -1) to (100 50 1)
1 by 1 by 1 MPI processor grid
reading atoms ...
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
reading bonds ...
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
read_data CPU = 0.004 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.