[lammps-users] couple lammps with quantum espresso

Hello all,
I follow the example from couple/lammps_quest. Instead of quest
package the pwscf
was used with corresponding modifications. The resulting code is working with
dynamics and dempfered dynamic. However I can not minimize the structure.

The output is:

LAMMPS (23 Sep 2010)
Lattice spacing in x,y,z = 5.39761 5.39761 5.39761
Created orthogonal box = (0 0 0) to (5.39761 5.39761 5.39761)
  1 by 1 by 1 processor grid
Created 8 atoms
WARNING: Resetting reneighboring criteria during minimization
ERROR: Fix external callback function not set

As an lammps input I used script:

units metal
dimension 3
atom_style atomic
atom_modify sort 0 0.0

lattice diamond 5.3976054
region box block 0 1 0 1 0 1
create_box 1 box
create_atoms 1 box
mass 1 28.08

velocity all create 10.0 87293 loop geom

fix 1 all external
dump 2 all custom 1 dump.md id type x y z fx fy fz

thermo 1

min_style cg

minimize 1.0e-6 1.0e-6 100 1000

Question: Does minimixation allows to use call_back fuction? And if
no, what shuld I change?

Thank you in advance
German Samolyuk

To use fix external with a minimization, a couple hooks
need to be added to the fix. I'll do it this AM.

But the error you are getting about the callback not being
set, is one your caller program has to take care of,
same as for dynamics.

Steve

German - I posted the patch for this - 22Nov10.
However, I realized something that may be problematic
about using this during a minimization. LAMMPS should
now get forces from the external program (quantum in this
caes), but will know nothing about energy associated
with those forces. So it may not work properly in a minimize.
I would try it with a 0.0 energy tolerance and see what
happens. It could be that there needs to be a 2nd callback
for energy or an extra arg to allow the external program
to also return an energy.

Steve

Steve,
Thank you.
Yes, I understand the problem with energy. Especially it we have two
regions, the forces from DFT for one region and
from MD for the second one.
German

Steve,
How I can run parallel DFT code "couple/lammps_quest".
If call system("mpirun -np 24 pw.x -npool 4 < pwscf.in >
pwscf.screen"), I slightly modified code to use quantum espresso,
will it execute it correctly? How I can run it if number of processors
should be initialized in batch file not in call of mpirun?
German

I'm not clear on what you are doing with espresso.
The LAMMPS/Quest example in couple has
a driver program lmpsqt.cpp which invokes LAMMPS
as a library and Quest as a stand-alone code (when
LAMMPS makes a callback to the driver). lmpqst
is what you would launch with mpirun or a batch
script. lmpsqt has a call to MPI_init() and queries
the number of procs it is running on via MPI_Comm_size.
It doesn't matter whether it is launched by mpirun or a batch script.

Steve

Steve,
I do exactly the same.
"lmpsqt.cpp which invokes LAMMPS
as a library and pwscf (qespresso code) as a stand-alone code (when
LAMMPS makes a callback to the driver)".
However in the LAMMPS/Quest example the lammps part runs parallel, the
Quest (pwscf) code not.
My question is how to run quest code in parallel also?
German

My question is how to run quest code in parallel also?

I assume you mean qespresso in parallel. Ideally the
quantum code would have an interface like LAMMPS.
I.e. a library interface that allows the driver program
to invoke it either in serial or parallel. Quest is not
(yet) set up that way, so that is why it is invoked
via a system() call as:
    system("lcao.x > lcao.screen");

If you wanted to do this invocation in parallel,
you might be able to do something like:

if (me == 0) system("mpirun -np 4 qespresso ...");
MPI_Barrier(MPI_COMM_WORLD);

But I'm not sure that would work. I.e. would the
other procs and the OS allow qespresso to run on them while
they wait on the barrier. Haven't tried it.

The right way to do this is for the target code
to have a library interface that is designed
to allow multiple processors to call it simultaneously, while
passing in the MPI communicator.

Steve

I asked an MPI expert here about this. He
said what I suggested below will probably
work on some parallel machines like your
multi-core workstation or that run full Linux
on every processor. But will not work on
many bigger parallel machines that have
service and compute nodes or that restrict
what the OS on a compute node can do. E.g.
the Cray XT series.

He tole me there is also, in MPI 2.0, a MPI_Spawn
command that formalizes what I suggested
below. I.e. it allows you to fork off a new MPI
job from a parent and gives you a communicator
to talk to the child job. But similar restrictions
apply. E.g. it would work on a desktop machine,
but not on something like the Cray XTs.

So the ultimate solution is what I said: have
any child code called by the driver be invoked
thru a library interface. That will work on any
machine.

Steve

Steve,
Yes. it looks like I should do it through the library interface.
Thanks
German

Steve,
In the code lammps_quest how I can access atoms positions and forces
of particular group only?
German

Steve,
In the code lammps_quest how I can access atoms positions and forces
of particular group only?

german,

you have a pointer to the lammps object that is
being created, that in turn has a pointer to the
group object from which you can find the group
id and via the group bitmask identify the atoms
that are in the group. just like it is done
elsewhere within lammps.

cheers,
    axel.

Axel,
Thanks!
This is how I'm trying to do it.

int groupbit = lmp->group->bitmask[0];

However the compiler return error

make -f Makefile.openmpi
mpic++ -g -O -DMPICH_IGNORE_CXX_SEEK -I../library -M lmppwscf.cpp > lmppwscf.d
mpic++ -g -O -DMPICH_IGNORE_CXX_SEEK -I../library -c lmppwscf.cpp
lmppwscf.cpp(97): error: pointer to incomplete class type is not allowed
    int groupbit = lmp->group->find("GDFT");
                   ^

compilation aborted for lmppwscf.cpp (code 2)

Don't see what I do wrong.

Best,
German

Axel,
Thanks!
This is how I'm trying to do it.

int groupbit = lmp->group->bitmask[0];

However the compiler return error

make -f Makefile.openmpi
mpic++ -g -O -DMPICH_IGNORE_CXX_SEEK -I../library -M lmppwscf.cpp > lmppwscf.d
mpic++ -g -O -DMPICH_IGNORE_CXX_SEEK -I../library -c lmppwscf.cpp
lmppwscf.cpp(97): error: pointer to incomplete class type is not allowed
int groupbit = lmp->group->find("GDFT");
^

compilation aborted for lmppwscf.cpp (code 2)

Don't see what I do wrong.

not including the proper header files with the full class definitions?

axel.

Axel,
It works. I should be more accurate.
Thanks
German

Dear German,

So am I wrong to assume from your last statement that you seem to have PWSCF working?

How soon until you submit it?

Will you make an advanced version available?

Nice work.

Best.

Josh

dear josh,

writing a proper qm/mm coupling requires much more than
making two codes talk to each other. that is only the beginning.
particularly plane wave based DFT codes have unexpected
quirks (you need to decouple the periodic images of the QM
system, you have to avoid "spill-out", you have to implement
a scheme to terminate "dangling bonds").

cheers,
     axel.

Dear Josh,

It's at the beginning stage. I'll let you know as soon as it will do
something reasonable.

Alex,
I have a question to you. The library function 'void Many2One::gather"
from the directory couple/library
modifies indexes of arrays. As a result
the cycle over atomic positions with condition

if (info->lmp->atom->mask[i] & groupbit) {

return wrong set of atomic position
How can this problem can be overcomed?
thank you in advance
German

Dear Josh,

It's at the beginning stage. I'll let you know as soon as it will do
something reasonable.

Alex,
I have a question to you. The library function 'void Many2One::gather"
from the directory couple/library
modifies indexes of arrays. As a result
the cycle over atomic positions with condition

if (info->lmp->atom->mask[i] & groupbit) {

return wrong set of atomic position

this code segment doesn't make any sense.
you cannot loop over distributed data and then
assume that it would work with data that has
been gathered from all nodes.

How can this problem can be overcomed?

you'd have to implement a version of Many2All::gather
for lists of integers where you'd combine and sort the
group mask data in the same way you collect positions.

i don't know what you want to use this for, but my gut feeling
is that you are trying to do something the wrong way around.

axel.

I'd like to determine group of atoms to work within DFT and rest of
atoms through MD, also boundary atoms. Hoe.
I definitely can simply check if atoms are inside specific region.
However for me it was more comfortable to work with
different groups of atoms.
German