Is there a way I could call for the total pair-force between two atoms for a compute regardless of the pair style used in the simulation?
Is there a way I could call for the total pair-force between two atoms for a
compute regardless of the pair style used in the simulation?
what exactly do you mean by "total pair force"?
and in what way does that have to be independent of the pair style in use?
axel.
By total pair force I mean the sum of all pair forces defined by the user between an atom i and an atom j. For example had I defined a hybrid use of LJ and coulomb I would want the coulomb force and LJ force between i and j either summed or separate so I can sum in the compute. As for "independent of pair style in use" I meant through the interfaces and data structures of LAMMPs so I don't have to code the compute calling specific pair styles.
By total pair force I mean the sum of all pair forces defined by the user between an atom i and an atom j. For example had I defined a hybrid use of LJ and coulomb I would want the coulomb force and LJ force between i and j either summed or separate so I can sum in the compute. As for "independent of pair style in use" I meant through the interfaces and data structures of LAMMPs so I don't have to code the compute calling specific pair styles.
if you want to compute pairwise force from the different component of
a hybrid pair style manually, you will have to do something along the
lines of the following:
1) add your compute class as friend class to pair_hybrid.h
2) add to the appropriate places the code segments below.
#include "pair_hybrid.h"
[...]
PairHybrid *hybrid = NULL;
if (strncmp(force->pair_style,"hybrid",6) == 0) {
hybrid = (PairHybrid *)force->pair;
} else error->all(FLERR,"This compute will only work with a hybrid pair style");
[...]
double epair, fpair;
for (int istyle = 0; istyle < hybrid->nstyles; ++istyle) {
epair = fpair = 0.0;
if (hybrid->styles[istyle]->single_enable) {
epair = hybrid->style[istyle]->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
}
[...]
}
now you compute in each loop iteration the individual contributions of
energy and force for a given i,j pair, for all hybrid sub-styles
providing the Pair::single() function.
axel.
Pair hybrid has a single() method you can call from a compute.
(this is what compute group/group does). That method
has 2 args for factor_coul and factor_lj which are prefactors
on the LJ and Coul terms. If you call it twice with
either factor = 0.0, you should get the LJ and Coul contributions
to force and energy.
PairHybrid does its own checking that you are using it’s single()
method with sub-styles that support it.
Steve
Thanks, the compute works as intended with that. Unfortunately the supported pair styles seems to be limited for the single() method. I’m using an overlay of buckingham and coulomb and its reporting a lack of support, is the only way out to implement a single method for every pair style I might use?
So far as I know every pair style in LAMMPS that is a true pair interaction (as opposed to many-body),
has a single() method defined. Both the Buckingham and cou/cut or coul/long support it.
Steve
looks like hybrid/overlay contains so single() definition, is that why?
It derives from PairHybrid which does the single settings, depending on the
sub-styles. Unless you post an input script using a command
that uses the single() method, like compute group/group, no one can
tell if the problem is with LAMMPS or your script.
Steve
Here is the input script, the call is made in the compute with id “srf” :
SrTiO3 Thermal simulation
units metal
dimension 3
boundary p p p
atom_style charge
neighbor 2.5 bin
neigh_modify delay 5
read_data data.bdry
read_restart restart.1280000
Potential
pair_style hybrid/overlay buck 9.0 coul/wolf 0.26506 9.0
pair_coeff 1 3 buck 1769.51 0.319894 0.0
pair_coeff 2 3 buck 14567.4 0.197584 0.0
pair_coeff 3 3 buck 6249.17 0.231472 0.0
pair_coeff * * coul/wolf
dielectric 1.0
define groups
region 1 block EDGE EDGE EDGE EDGE -20.0 20.0 units box
region 2 block EDGE EDGE EDGE EDGE 791.6 EDGE units box
region 3 block EDGE EDGE EDGE EDGE EDGE -791.6 units box
group source region 1
group sink region 2
group sink2 region 3
initial velocities
compute 1 all temp
compute 2 source temp
compute 3 sink temp
compute KE all ke/atom
compute srf all surfstress z 200
fix 1 all nvt temp 600 600 0.055
fix 3 all temp/rescale 100 0.01 0.01 0.002 1.0
variable kB equal 8.6171092e-5
variable t equal f_7[1][3]/1.5/v_k
variable k equal f_7[1][3]/1.5/c_1
variable t atom c_4/1.5/v_k
fix 7 all ave/spatial 70 1 70 z lower 0.1 v_t units reduced file tp.txt
fix 9 all ave/spatial 1 4000 4000 z lower 4.0 v_t units box file tpwindow1.txt
fix 10 all ave/spatial 1 4000 4000 z lower 4.0 v_t units box file tprunning1.txt ave running
timestep 0.0002
thermo 100
thermo_style custom step temp c_srf
run 100000
So what command is giving an error about single() not being suppported
by the pair style?
Steve
compute surfstress calls single().
I don’t want to debug a compute you wrote. Please post a simple
as possible input script that doesn’t work with compute group/group
due to issues with pair hybrid and single(). If you can’t make
that fail, then the problem is in your new compute …
Steve
looks like pair_style coul/wolf has Pair::single_enable set to 0.
how about trying coul/cut or coul/long?
axel.
It works with compute group/group as Steve suggested which also calls single(), there’s some subtle issue here. I practically xerox copied compute group/group to create the compute I made and just deleted all the group2 arg data members and assignments in the implementation and made no modification to single() so I’m not sure whats causing it yet.
It works with compute group/group as Steve suggested which also calls
single(), there's some subtle issue here. I practically xerox copied compute
group/group to create the compute I made and just deleted all the group2 arg
data members and assignments in the implementation and made no modification
to single() so I'm not sure whats causing it yet.
this would only make sense, if you are using a LAMMPS version older
than 1 Feb 2014, which is the patch level where the single setting in
coul/wolf was disabled.
axel.
The version I modified is Feb 2015.
then you may have defined compute group/group but are not actually
using it, e.g. access its data.
otherwise you should get:
Setting up Verlet run ...
Unit style : real
Current step : 0
Time step : 2
Memory usage per processor = 7.29551 Mbytes
Step Temp PotEng KinEng Press 1
ERROR on proc 0: Pair hybrid sub-style does not support single call
(../pair_hybrid.cpp:783)
ERROR on proc 2: Pair hybrid sub-style does not support single call
(../pair_hybrid.cpp:783)
ERROR on proc 3: Pair hybrid sub-style does not support single call
(../pair_hybrid.cpp:783)
ERROR on proc 1: Pair hybrid sub-style does not support single call
(../pair_hybrid.cpp:783)
Yea, you’re right I tried group all on all and got the same error, although I was accessing its data before. I had specified two arbitrary groups and the data was just zero since the groups were so far apart they had no neighbor interaction. Thanks for the help.
That’s just an error. Someone probably copy/pasted to create the coulwolf
sytle and didn’t see the setting. Coul/wolf provides a single method, so just
delete single_enable = 0 from the constructor.
Steve