Fix adapt keyword "scale" is not implemented fully

Dear everyone,

I found that in fix adapt the keyword "scale" is not functioning for all
the attributes that are offered.

Let's say I want to change the charges of the atoms and scale them by a
factor. Normally one would use "scale yes". But I found that it does not
matter whether I set "scale yes" or "scale no" - it will always set the
charges to the value I input (behaving like „scale no“), and never just
multiply them with this value.

So I looked into fix_adapt.cpp, more precisely into „change_settings()“,
and found that the scaleflag is only used for the attribute „pair“. But
for „atom“ (and also „kspace“) there is no if-clause handling the
scaleflag. Let‘s take the example of „atom charge“, which I need the
most currently:

   } else if (ad->which == ATOM) {

     // …

     if (ad->aparam == DIAMETER) {

    // …

     } else if (ad->aparam == CHARGE) {
       double *q = atom->q;
       int *mask = atom->mask;
       int nlocal = atom->nlocal;
       int nall = nlocal + atom->nghost;

       for (i = 0; i < nall; i++)
         if (mask[i] & groupbit) q[i] = value;
     }
   }

So I thought of implementing it myself. Of course it is easy to just
write a simple if-clause like this:

       if (scaleflag) {
          for (i = 0; i < nall; i++)
            if (mask[i] & groupbit) q[i] = // insert calculation for
scaled charge
       } else {
         for (i = 0; i < nall; i++)
           if (mask[i] & groupbit) q[i] = value;
       }

But I have a problem with the calculation: I cannot just write q[i] =
q[i]*value because then it would multiply the charges iteratively
everytime when the fix is called. Lets say value is a number between 0
and 1 then the charges would get smaller and smaller...
I have to code something like q[i] = q_start_of_run[i]*value; but my
problem is how to get q_start_of_run[i], which is the value at the
beginning of the simulation run. I think the original charges are stored
in fix_chg→vstore (see restore_settings() and post_constructor() ) so I
need to base my calculation on this, but all my tries so far have
resulted in memory access errors.
How to do this properly?

Finally, I want to remark that these problems also apply to fix
adapt/fep.
Furthermore, in the Lammps version 31Mar2017, I cannot find where the
keyword „bond“ is ever used in fix adapt. This function is suggested in
the online manual http://lammps.sandia.gov/doc/fix_adapt.html but the
code of fix_adapt.cpp does not contain a „bond“ option.

Thank you for your consideration.

Best wishes,

Markus Koch

Dear everyone,

I found that in fix adapt the keyword "scale" is not functioning for all
the attributes that are offered.

Let's say I want to change the charges of the atoms and scale them by a
factor. Normally one would use "scale yes". But I found that it does not
matter whether I set "scale yes" or "scale no" - it will always set the
charges to the value I input (behaving like „scale no“), and never just
multiply them with this value.

So I looked into fix_adapt.cpp, more precisely into „change_settings()“,
and found that the scaleflag is only used for the attribute „pair“. But
for „atom“ (and also „kspace“) there is no if-clause handling the
scaleflag. Let‘s take the example of „atom charge“, which I need the
most currently:

   } else if (ad->which == ATOM) {

     // …

     if (ad->aparam == DIAMETER) {

    // …

     } else if (ad->aparam == CHARGE) {
       double *q = atom->q;
       int *mask = atom->mask;
       int nlocal = atom->nlocal;
       int nall = nlocal + atom->nghost;

       for (i = 0; i < nall; i++)
         if (mask[i] & groupbit) q[i] = value;
     }
   }

So I thought of implementing it myself. Of course it is easy to just
write a simple if-clause like this:

       if (scaleflag) {
          for (i = 0; i < nall; i++)
            if (mask[i] & groupbit) q[i] = // insert calculation for
scaled charge
       } else {
         for (i = 0; i < nall; i++)
           if (mask[i] & groupbit) q[i] = value;
       }

But I have a problem with the calculation: I cannot just write q[i] =
q[i]*value because then it would multiply the charges iteratively
everytime when the fix is called. Lets say value is a number between 0
and 1 then the charges would get smaller and smaller...
I have to code something like q[i] = q_start_of_run[i]*value; but my
problem is how to get q_start_of_run[i], which is the value at the
beginning of the simulation run. I think the original charges are stored
in fix_chg→vstore (see restore_settings() and post_constructor() ) so I
need to base my calculation on this, but all my tries so far have
resulted in memory access errors.
How to do this properly?

​first of all, scaling of coulomb interactions can be done by adjusting the
"scale" property of coulomb pair styles. no need to explicitly scale
charges. ​

regardless of that, ​what you want to implement does not require any
changes to the source code. but you should use the "scale no" option and
instead define an atom style variable that computes the current value you
want to be used based on the time step number. the original charge
information can be stored via fix store/state and then referenced in your
atom style variable.
and then also used to restore the original value with set.

Finally, I want to remark that these problems also apply to fix
adapt/fep.
Furthermore, in the Lammps version 31Mar2017, I cannot find where the
keyword „bond“ is ever used in fix adapt. This function is suggested in
the online manual http://lammps.sandia.gov/doc/fix_adapt.html but the
code of fix_adapt.cpp does not contain a „bond“ option.

​please keep in mind, that the online documentation always reflects the
very latest LAMMPS version (13 April 2017 as of today) and then have a look
at: ​http://lammps.sandia.gov/bug.html

​axel.​

Dear Axel,

thank you for your fast reply.

One option I have explored but it did not work for me:

  1. Actually I would also prefer to scale only the coulomb interactions, but for that I encounter a different problem.
    Let me explain shortly what happens. I simulate a molecule in water using the class2 force field PCFF. I normally use pair_style lj/class2/coul/cut for my simulation but I want to scale the LJ and the Coulomb interactions independently (first turn on LJ from none to fully on, then the Coulomb from none to full interactions).

I have already implemented the extract() method in all the pair_lj_class2 styles and this works fine. To apply fix_adapt independently to LJ and Coul. interactions I figured (as was also suggested in some threads here) one could use “pair_style hybrid/overlay” with lj/class2 and coul/cut and then reference them independently in the fix adapt. This works fine, but a problem occurs when I want to use compute pe/tally to evaluate the coulomb and LJ energies independently (for the water-water, molecule-molecule and molecule-water interactions).

Using this hybrid/overlay pair style and then trying to output the energies from pe/tally, I get an error “Energy was not tallied on needed timestep”.
This did not occur when I was using lj/class2/coul/cut with pe/tally. All I changed between scripts were the lines regarding pair styles and pair coefficients. Otherwise, the simulation with hybrid/overlay runs fine (i.e. if I don’t compute and output data from pe/tally).

Is it expected that hybrid/overlay and pe/tally cannot be used at the same time? Or is there some other mistake I am doing (my input script is poste at the end of my message)?

I want to avoid using rerun to get the energy contributions, which I tested before (it works) but I’d prefer to produce the data on the fly.

  1. The option to scale the charges with fix store/state and fix adapt with “scale no” I will test next. Thank you.

Thank you also for reminding me that the manual is for the newest version, not the newest stable one. I didn’t notice that before but now I see that it is also noted in the header of manual pages.

Here is my input script for the problem when using pair_style hybrid/overlay and compute pe/tally together, maybe its something in there?

units real
dimension 3
atom_style full
boundary p p p
special_bonds lj/coul 0 0 1
neigh_modify one 5000

pair_style hybrid/overlay lj/class2 20.0 coul/cut 20.0
###pair_style lj/class2/coul/cut 20.0 20.0

bond_style class2
angle_style class2
dihedral_style class2
improper_style class2

read_data minimized.data

I wrote all pair_coeffs explicitly just to be sure

everything is set properly, all the rest is defined in

minimized.data

pair_coeff 1 1 lj/class2 0.064 4.01 20
pair_coeff 1 2 lj/class2 0.0627813 4.06739 20
pair_coeff 1 3 lj/class2 0.0644341 4.04056 20

pair_coeff 10 11 lj/class2 0.00336155 3.21479 20
pair_coeff 11 11 lj/class2 0.274 3.608 20
pair_coeff 1 1 coul/cut
pair_coeff 1 2 coul/cut

pair_coeff 10 11 coul/cut
pair_coeff 11 11 coul/cut

group star type 1 2 3 4 5 6 7 8 9
group water type 10 11

timestep 1.0 #fs

velocity all create 300.0 921841 rot yes dist gaussian
fix 1 all nvt temp 300.0 300.0 100.0

just some example values

variable lj equal 1.0
variable ch equal 0.5

fix adapt_lj all adapt 1 pair lj/class2 epsilon 19 1011 v_lj scale yes reset no
fix adapt_coul all adapt 1 pair coul/cut scale 19 1011 v_ch scale yes reset no

compute e_all all pe/tally all
compute e_star star pe/tally star
compute e_water water pe/tally water
compute e_int star pe/tally water

compute red_all all reduce sum c_e_all[1] c_e_all[2]
compute red_star all reduce sum c_e_star[1] c_e_star[2]
compute red_water all reduce sum c_e_water[1] c_e_water[2]
compute red_int all reduce sum c_e_int[1] c_e_int[2]

output interaction en. of star and water - 1 is VdW, 2 is Coul

thermo_style custom step c_red_int[1] c_red_int[2]

thermo 5

run 10

Dear Axel,

thank you for your fast reply.

One option I have explored but it did not work for me:

1) Actually I would also prefer to scale only the coulomb interactions,
but for that I encounter a different problem.
Let me explain shortly what happens. I simulate a molecule in water using
the class2 force field PCFF. I normally use pair_style lj/class2/coul/cut
for my simulation but I want to scale the LJ and the Coulomb interactions
independently (first turn on LJ from none to fully on, then the Coulomb
from none to full interactions).

I have already implemented the extract() method in all the pair_lj_class2
styles and this works fine. To apply fix_adapt independently to LJ and
Coul. interactions I figured (as was also suggested in some threads here)
one could use "pair_style hybrid/overlay" with lj/class2 and coul/cut and
then reference them independently in the fix adapt. This works fine, but a
problem occurs when I want to use compute pe/tally to evaluate the coulomb
and LJ energies independently (for the water-water, molecule-molecule and
molecule-water interactions).

Using this hybrid/overlay pair style and then trying to output the
energies from pe/tally, I get an error "Energy was not tallied on needed
timestep".
This did not occur when I was using lj/class2/coul/cut with pe/tally. All
I changed between scripts were the lines regarding pair styles and pair
coefficients. Otherwise, the simulation with hybrid/overlay runs fine (i.e.
if I don't compute and output data from pe/tally).

Is it expected that hybrid/overlay and pe/tally cannot be used at the same
time? Or is there some other mistake I am doing (my input script is poste
at the end of my message)?

​i will have to set up a test system for myself and track this one down. i
don't think i have ever tested the code in USER-TALLY with hybrid/overlay.
i can see, why there might be a conflict.

I want to avoid using rerun to get the energy contributions, which I
tested before (it works) but I'd prefer to produce the data on the fly.

2) The option to scale the charges with fix store/state and fix adapt with
"scale no" I will test next. Thank you.

​please keep in mind, that when scaling the charges, you will apply
multiply with the scale factor twice for each pair compared to scaling the
coulomb pair interaction and you'll have to account for that (e.g. by using
sqrt()), if you want to get comparable results.

​axel.​

​[...]

Using this hybrid/overlay pair style and then trying to output the

energies from pe/tally, I get an error "Energy was not tallied on needed
timestep".
This did not occur when I was using lj/class2/coul/cut with pe/tally. All
I changed between scripts were the lines regarding pair styles and pair
coefficients. Otherwise, the simulation with hybrid/overlay runs fine (i.e.
if I don't compute and output data from pe/tally).

Is it expected that hybrid/overlay and pe/tally cannot be used at the
same time? Or is there some other mistake I am doing (my input script is
poste at the end of my message)?

​i will have to set up a test system for myself and track this one down. i
don't think i have ever tested the code in USER-TALLY with hybrid/overlay.
i can see, why there might be a conflict.

​indeed, some small changes to the sources are needed to ​properly
propagate the callback registration to the hybrid substyles.

​i am attaching a patch to the latest sources with the required changes.
they are also in a pending pull request and will likely be included into
the next (unstable​) release version and the git master branch.

axel.

add-user-tally-for-hybrid.diff.gz (742 Bytes)

Thank you very much, looking forward to the patch!

Thank you very much, looking forward to the patch!

if desired, ​you can download a snapshot from my pull request branch here: ​
https://github.com/akohlmey/lammps/tree/small-fixes-and-updates

​axel.​