How to use two different molecular templates in fix-gcmc

LAMMPS LAMMPS Mailing List Mirror
1

Dear lammps users,

I want to do gcmc where two kind of molecules are exchanged, say, H2 and CH4.

In the current lammps code, it is not allowed to do this, thus I have to do some modifications.

I plan is to add a new class, fix_gcmc_plus, then in the script, I will have

fix 1 H2 gcmc 10 10 10 …mol h2 …

fix 2 CH4 gcmc_plus 10 10 10 …mol ch4 …

The problem is, I guess, only one molecular template is allowed in lammps code. The molecular information for which molecule to be inserted, is saved in a class named “onemols”. It seems, for both fix_gcmc and fix_gcmc_plus, they want to read the “onemol”. As a result, I always got a segmentaltion fault.

Can somebody give me some guides where and how to modify the code? Thanks in advance!

Best,

Yongbiao

2

Dear lammps users,

    I want to do gcmc where two kind of molecules are exchanged, say, H2 and
CH4.
    In the current lammps code, it is not allowed to do this, thus I have to
do some modifications.
    I plan is to add a new class, fix_gcmc_plus, then in the script, I will
have

that is a bad plan. the fix gcmc documentation states you can already
can exchange multiple molecules if you use multiple gcmc fixes.

    fix 1 H2 gcmc 10 10 10 .......mol h2 ......
    fix 2 CH4 gcmc_plus 10 10 10 .....mol ch4 ......

   The problem is, I guess, only one molecular template is allowed in lammps
code. The molecular information for which molecule to be inserted, is saved
in a class named "onemols". It seems, for both fix_gcmc and fix_gcmc_plus,
they want to read the "onemol". As a result, I always got a segmentaltion
fault.

onemols is just a pointer to atom->molecule[] which is the list of
possible molecule templates added with the "molecule" command.
have you actually tried to understand where the segfault is coming
from? without that skill, you will have a very hard time do
modifications at a rather complex piece of clode.

the reason for the segfault is an inconsistency (one might call it a
bug) in how "onemols" is initialized and then later used. a careful
look at the source code and my change reveals, that the it will
*always* cause segfaults, when a different molecule template than the
very first one defined is used. regardless whether you have one or
more fix gcmc instances defined.
the following change corrects that.

diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index a9c7ad4..b0c1982 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -253,8 +253,8 @@ void FixGCMC::options(int narg, char **arg)
         error->warning(FLERR,"Molecule template for "
                        "fix gcmc has multiple molecules");
       mode = MOLECULE;
- onemols = &atom->molecules[imol];
- nmol = onemols[0]->nset;
+ onemols = atom->molecules;
+ nmol = onemols[imol]->nset;
       iarg += 2;
     } else if (strcmp(arg[iarg],"region") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");

  Can somebody give me some guides where and how to modify the code? Thanks
in advance!

you need to spend more time studying C++ programming (i assume you
already have a suitable book that you learned C++ from) and practice
reading, understanding, and debugging C++ code.
here is a URL for a book, that is particularly useful for practicing
reading and understanding code (most regular programming books don't
do that, they mostly only teach syntax).

http://www.amazon.com/How-Not-Program-Programs-Working/dp/1886411956/ref=sr_1_1?ie=UTF8&s=books&qid=1237162335&sr=1-1

debugging is like solving integrals: it needs a lot of practice.

axel.

3

Thanks Axel. This was indeed an error in the code that had no effect
in the most common use case where the molecule template specified by
fix gcmc was also the only one appearing in the input script. Anyone
wanting to insert two different molecules using two different fix gcmc
commands will need to update fix_gcmc.cpp with Axel's fix.

4

Thanks Axel. Your improvement really helps.

But I am afraid that more modification should be done.

I made tests after Axel’s fix, and found that it can run several thousands of steps, then

a new “segmentation faul” came out without more information.

I am trying to debug it, which may take much time because I am a beginner at c++.

If somebody can do tests, too, I will be very appreciate!

Thanks all of you!

Yongbiao

5

Please post a small-as-possible example input file demonstrating the problem.

2015-11-11 18:44 GMT-07:00 Yongbiao Yang <[email protected]...>:

6

Thanks for you kind help.

The mainly error information is as follows:
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (…/neighbor.cpp:441)
WARNING: System is not charge neutral, net charge = 1.6952 (…/kspace.cpp:297)
ERROR: Lost atoms: original 2850 current 2848 (…/thermo.cpp:395)

For details, see the input files as follows (the input data file for water is in the attachment)

----------------------------------------------------------input--------------------------------------------------------------------

units real
atom_style full
boundary p p p
pair_style lj/cut/coul/long 10.0
pair_modify tail yes
bond_style harmonic
angle_style harmonic
special_bonds lj/coul 0 0 0.5
kspace_style pppm 0.00001
variable molname index yyb
variable equil index heat
variable itemp index 300
variable ftemp index 300

read_restart ddrst.rst2

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

group water type 1 2
molecule hho …/00model/hho.txt offset 0 0 0 0 0

info group

log log.{molname} ##_{ensemble}_${nanosec}

timestep 1.0

#fix fixshake all shake 0.0001 20 0 b 1 a 1
fix fixnvt all nvt temp {itemp} {ftemp} 50.0
fix fixmomentum all momentum 1 linear 1 1 1
fix fixrecenter all recenter INIT INIT INIT

fix fixgcmc water gcmc 2 10 10 0 777 ${ftemp} 0.0 0.1 mol hho maxangle 90 full_energy

thermo 2
thermo_style custom step time atoms temp press lx density

restart 50 ddrst.rst1 ddrst.rst2

dump dumptrj all atom 100 {molname}_all.trj ##{ensemble}_${nanosec}.trj
dump_modify dumptrj image yes scale yes

run 10000 ## 1000000 # for upto ${nanosec}

unfix fixnvt # unfix fixmomentum
unfix fixrecenter
unfix fixgcmc
undump dumptrj

----------------------------------------------------------input--------------------------------------------------------------------

----------------------------------------------------------output--------------------------------------------------------------------

TACC: Setting memory limits for job 2903494 to unlimited KB
TACC: Dumping job script:

in.data (569 KB)

in.parameter (1.33 KB)

log.out (6.29 KB)

7

ERROR: Lost atoms: original 2850 current 2848 (../thermo.cpp:395)

This is not a "segmentation fault" and it does not indicate a bug in
the code. We call this a LAMMPS error message. Probably you have
inserted a molecule in a high energy (and large force) configuration.
As a result, two atoms moved a large distance, and became "lost." If
you are operating at a thermodynamically reasonable chemical
potential, the probability of this happening is much less than
1/10,000.

Aidan

2015-11-12 14:32 GMT-07:00 Yongbiao Yang <[email protected]...>: