Problems with fix addforce and referring to variables

Hi Everyone,

I am using the (21 Dec 2011) LAMMPS version and am trying to apply a potential dependent on only the x coordinate of particles in my system using the fix addforce command, but have run into a few problems. After searching the mail list and the user manual for what I believe were relevant keywords and commands, I am still rather confused about what I should do. I apologize for the length of this message, but I tried to keep it as simple and concise as possible.

To verify that the potential has the shape I want, I decided to:

  1. create a box and insert a single test particle
  2. loop over several x coordinates for that particle. For each iteration, I use “run 0” then print the x-coord and the potential energy to an output file.
  3. plot the potential energy vs. x-coord
    The following is an excerpt of my input script:

Comments below.

Steve

Hi Everyone,

I am using the (21 Dec 2011) LAMMPS version and am trying to apply a
potential dependent on only the x coordinate of particles in my system using
the fix addforce command, but have run into a few problems. After searching
the mail list and the user manual for what I believe were relevant keywords
and commands, I am still rather confused about what I should do. I apologize
for the length of this message, but I tried to keep it as simple and concise
as possible.

To verify that the potential has the shape I want, I decided to:

create a box and insert a single test particle
loop over several x coordinates for that particle. For each iteration, I use
"run 0" then print the x-coord and the potential energy to an output file.
plot the potential energy vs. x-coord

The following is an excerpt of my input script:
____________________

variable force_x atom 2x
fix the_potential sodium addforce v_force_x NULL NULL
fix_modify the_potential energy yes
variable curr_x equal 5
variable delta_x equal 1
variable valx loop 5000
        label next_one
        variable curr_x equal \{curr\_x\}\-{delta_x}
        set group sodium x \{curr\_x\} &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;run 0 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;variable spe equal f\_the\_potential &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;shell echo {curr_x} \{spe\} &gt;&gt; pe\_vs\_x\.dat &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if &quot;{curr_x} <= 0" then "jump in.no_crystal_edit
exit_loop"
        next valx
        jump in.no_crystal_edit next_one
label exit_loop

____________________

I receive the following error message in my log file after line 10: "ERROR:
Invalid syntax in variable formula"

My questions:

* According to my version's LAMMPS manual fix addforce doc page, it appears
I need to use v_variable, not ${variable} to call my force variables.
Additionally, v_variable must be of the "atom" type if I want to have a
per-atom potential calculated in the future. Am I correct in my
interpretation?

yes

* Why is the following the case?

- When I hard-code "2x" into line 2 and change nothing else in the script,
my output file (pe_vs_x.dat) file gives me exactly what I would expect:

5 -10
4 -8
3 -6
2 -4
1 -2
0 0

- But if I hard-code something like "sin(x)" or "exp(x)," the potential
energy for every x-coord is 0 in my output file. I also tried "sin(2.0*x)"
to make the numbers floats in case that was the problem, but still got 0 PE
for all x-coord.

* Do I need to modify line 11 of the excerpt in some way?

* Do you have any other suggestions for where I could look or what I could
do?

Thank you so much for your time and attention. I really appreciate any
thoughts or suggestions you have to offer!

With best regards,
Carol (incoming first-year graduate student at UCSB)

2x is in not a valid syntax for a variable formula. It has to be 2*x.
When I put these lines in bench/in.lj

variable force_x atom 2x
fix the_potential all addforce v_force_x NULL NULL
dump 1 all custom 10 tmp.dump id x y x fx fy fz
run 0

I get this error:
ERROR: Invalid syntax in variable formula (variable.cpp:713)

If you don't get that error, are you running the current version of LAMMPS?

If I replace 2x by 2*x, it runs fine and the dump file shows
that every atom has a force that is twice its x coord, which is
the correct result. I suggest you dump
forces on your atom into a file and look at them directly.
Once you are confident those are correct you can debug the
rest of your script and its loop to see if you can get them printed
to the screen.

Steve

Hi Steve,

Thank you for the suggestions and help! The temporary dump file you suggested is giving the forces I expected. However, when I try to include the potential energy, I get an error:

1. variable force_x atom 2*x
2. fix the_potential all addforce v_force_x NULL NULL
3. fix_modify the_potential energy yes
4. compute new_pe all pe
5. dump 1 all custom 10 tmp.dump id x y x fx fy fz c_new_pe
ERROR: Dump custom compute does not compute per-atom info (dump_custom.cpp:1149)

I can't find this error anywhere in the LAMMPS errors section. According to my version of the LAMMPS manual (Dec 2011) the dump custom command /should/ handle per-atom info. I believe that line 4 has the same syntax as that on the compute command manual page and that line 5 refers to my compute command appropriately. Am I mistaken? Or am I trying to obtain the new potential energy in an inappropriate way?

Best regards,
Carol

The error is listed in doc/Section_errors.html, like all error messages.

You need to use compute pe/atom (for anything you want to dump).
Compute pe calculates the global pot eng of the entire system.
The doc pages for the 2 computes explain this in some detail, as
does the howto page they reference.

Steve