Assign compute angle/local values to a variable

Hi all,

I would like to run a simulation with bonded and angle interactions. In this simulation, I want to remove angle interactions once a critical angle limit has been reached, similar to bonds in fix bond/break.

In my simulations, I output all the angles and angle types through “compute angle/ local” and “compute property/local” commands. Then I check all the angles in the list whether they meet the threshold, and if any of them do, I change their angle coefficients (harmonic in my case) to 0. Here is the input script and data file for a simple case of 3 atoms :


building simulation input

units si
boundary p p p
atom_style angle

bond_style harmonic
angle_style harmonic

read_data debug_angle_data_file.txt

group support id 1 3
group load id 2

pair_style none
atom_modify sort 10 1
neighbor 1.000000 bin
neigh_modify every 10 check yes

timestep 0.000005

fix 1 load addforce 0 0 -100
fix 2 support setforce 0 0 0
fix 3 all nve

compute 1 all angle/local theta eng
compute 2 all property/local atype
compute PotEn_Angle all pe angle

thermo 1000
thermo_style custom step pe ke etotal temp
dump 1 all xyz 1000
dump 2 all local 1000 THETA.txt c_1[1] c_1[2] c_2

variable ang equal c_1[1]
variable ang_type equal c_2

label loopa
variable a loop 1000
print “A = $a”
run 1000

if “{ang} < 179" then & "angle_coeff {ang_type} 0 0”

next a
jump debug_angle.txt loopa


%building input file
3 atoms
2 bonds
1 angles

1 atom types
2 bond types
1 angle types

-50 50 xlo xhi
-50 50 ylo yhi
-50 50 zlo zhi


1 0.001

Bond Coeffs

1 1000000 0.1
2 1000000 0.1

Angle Coeffs

1 0.9 180


1 1 1 0 0 0
2 2 1 0.1 0 0
3 3 1 0.2 0 0


1 1 1 2
2 2 2 3


1 1 1 2 3

My problem is that I cannot assign values of a local compute to a variable (eg variable ang equal c_1[1], where compute 1 outputs local data), hence I cannot access the results of angular computes. I would really appreciate any suggestions or workarounds!

Konstantinos Keremidis
PhD Candidate
Concrete Sustainability Hub
Department of Civil and Environmental Engineering
Massachusetts Institute of Technology

Generally, this sort of thing is better done internally, i.e. by a ‘fix.’ Ideally, something this simple would be accomplished for example by a modification to ‘fix adapt.’

Barring that, although fix bond/react was designed for considerably more complex situations, it might be able to do what you want by defining the end atoms of the angle as the ‘bonding atoms,’ and triggering a ‘reaction’ once they exceed a certain distance.

If the angle returns below the critical angle limit, do you want the angular forces to be switched back on again?
Although this sounds like this isn’t what you want, I mention it because it can be be accomplished easily with “angle_style table”. (Just create a text file with a function that resembles a parabola that has been clipped.)