Water example to test the PIMD+QE setup/installation?

[MQS-mark]

Hi Yuki Nagai, Would it be possible to add a single and 4 molecules water example to this repository to test if the setup of PIMD with QE works?

I tried running PIMD+QE but getting this error: Error termination at: subroutine setatomchanges_MPI

I can see in the source code that the input.dat file is opened and the keyword is needed, but in the manual of PIMD there is no explanation what the keyword means and how it should be defined:

source/setpiparams_MPI.F90 (from line 199 on):

!***********************************************************************
      subroutine setatomchanges_MPI
!***********************************************************************
!=======================================================================
!
!     set up MD parameters
!
!=======================================================================
!-----------------------------------------------------------------------
!     /*   shared variables                                           */
!-----------------------------------------------------------------------

      use common_variables, only : &
     &   iounit, atom_change, natom, myrank_world

!-----------------------------------------------------------------------
!     /*   local variables                                            */
!-----------------------------------------------------------------------

      implicit none

      integer :: i, j, k, nline, ierr

      character(len=9) :: char

!-----------------------------------------------------------------------
!     /*    read mass scaling parameters                              */
!-----------------------------------------------------------------------

      if ( myrank_world .eq. 0 ) then

         atom_change(1:natom) = 'NONE'

         open ( iounit, file = 'input.dat' )

         call search_tag ( '<atom_change>', 13, iounit, ierr )

         read ( iounit, *, iostat=ierr ) char

         read ( char, *, iostat=ierr ) nline

         if ( ierr .eq. 0 ) then

            do k = 1, nline

               read ( iounit, *, iostat=ierr ) char, i, j

               if ( ierr .ne. 0 ) then

                  backspace( iounit )

                  read ( iounit, *, iostat=ierr ) char, i

                  j = i

               end if

               if ( ( i .le. 0 ) .or. ( i .gt. natom ) ) i = 1
               if ( ( j .le. 0 ) .or. ( j .gt. natom ) ) j = natom

               atom_change(i:j) = char

            end do

         else

            close( iounit )

            open ( iounit, file = 'input_default.dat' )

            call search_tag ( '<atom_change>', 13, iounit, ierr )

            read ( iounit, *, iostat=ierr ) char

         end if

         close( iounit )

      end if

      call my_mpi_bcast_int_0_world( ierr )

      call error_handling_MPI &
     &   ( ierr, 'subroutine setatomchanges_MPI', 29 )

      call my_mpi_bcast_char_1_world( atom_change, len(atom_change), natom )

      return
      end

I used this example from the PIMD webpage as a basis: https://ccse.jaea.go.jp/software/PIMD/pi.en.html

And this is how I modified the input.dat:

<method>
PIMD

<ensemble>
NVT

<bath_type>
MNHC

<temperature>
300.d0

<time_bath>
10.d0

<nnhc>
4

<dt>
0.25d0

<nstep>
10000

<corrections>
0  0

<nbead>
1

<ipotential>
QE

<np_beads>
-1

<qe_input_file_name>
4_h2o.in

<qe_output>
1 10

<iboundary>
0

<iread_exit>
10

<iprint_std>
1

<iprint_bond>
10

<iprint_eavg>
10

<iprint_rest>
10

<iprint_trj>
10

<iprint_xyz>
1

<iprint_dip>
100

<iprint_rdf>
10

The 4_h2o.in file is declared as follows:

&control
    calculation = 'scf'
    restart_mode='from_scratch'
    prefix = '4_h2o'
    outdir='./outdir'
    pseudo_dir = '../../pimd.2.6.1/lib/qe/pseudo'
/
&system
    ibrav=1, celldm(1) =9.86, nat=12, ntyp=2,
    ecutwfc=25, nosym=.true.
/
&electrons
    conv_thr=1e-8
/
&ions
/
ATOMIC_SPECIES
 H  1.00784  H.pbe-rrkjus.UPF
 O  15.999   O.pbe-rrkjus.UPF
ATOMIC_POSITIONS (angstrom)
O 1.6759513692647068 -0.6002853757772342 -1.6309012548029673 
H  1.0279586472537519 -0.9024715781296282 -2.26876660715738 
H  1.6152391970432078 -1.2290779155851816 -0.9104727361895202
O -1.5834963228239756 2.008349069216274 1.232856299827159
H  -1.5747845721704647 1.0714599470204051 1.4334707422714463 
H  -0.8617765374263848 2.128261335924955 0.6141278649996753 
O -1.3650917294404066 -0.565485109588074 2.1821923460132715
H  -1.980150002454046 -0.467032993221947 -2.3069381847013353
H  -1.350079743262641 -1.5072050335791736 2.00606312800933171
O 0.8016958070392737 2.1437494450214936 -0.1479543311321656
H  0.6999897857190609 1.6465013944331024 -0.9606557140636844
H  1.6273193428536925 1.830277846952749 0.22375067383692304 
K_POINTS (automatic)
  1 1 1 0 0 0

Would be great to manage to make some water examples for PIMD+QE together.

[cometscome]

Dear Dr. Mark Nicholas Jones,

Thank you for giving me a question. Could you tell me a version of PIMD? I will ask Dr. Shiga, the developer of PIMD. And I write an answer when I get it.

In example directory, there is a "water" directory, which might be useful for you.

Best, Yuki

[MQS-mark]

Hi Yuki, I'm using the most recent PIMD version 2.6.1.

I also wrote already to Motoyuki Shiga and Bo Thomsen via mail. I will include you in that email thread if the discussion develops further via mail. They also have the link to this issue page here so that there is a possibility of documenting my progress/experience with PIMD+QE.

Thank you for the info about the example folder. I was not aware of that. Will test the examples today and update this thread when I have worked things out or have some information to document/report.

Best, mark

[cometscome]

Hi Mark,

I asked Dr. Shiga. He said, please use

0 , since this tag is used only for changing atomic mass to set arbitrary value. And I think that you did not use the latest input_default.dat, which can be found in example dir. Please copy input_default.dat to your directory. Then, you can use default tags including atom_change. You can see what tags are in input_default.dat. p.s. I found that you wrote new comments just now. Thank you for telling me the version info. >I also wrote already to Motoyuki Shiga and Bo Thomsen via mail. I will include you in that email thread if the discussion >develops further via mail. They also have the link to this issue page here so that there is a possibility of documenting my >progress/experience with PIMD+QE. I see. I think it is better to add more examples with QE so I will add water case wit QE. Thank you for giving us good comments and suggestion. Best, Yuki

[MQS-mark]

Great! Will try this out. Thank you for all the help and this very useful repository!

[MQS-mark]

So I have now updated my project folder for a 4 water molecules example with the input_default.dat file from the examples folder. This solved the issue with , so the example on the webpage (https://ccse.jaea.go.jp/software/PIMD/pi.en.html) is outdated with respect to the input_default.dat file and the current version of pimd (2.6.1).

I run the following command:

mpirun -np 8 pimd.mpi.x

Now I receive the following error which is hard for me to decipher where the problem exactly is (I receive the same error when I try to run the Si/qe_pimd example):

*** An error occurred in MPI_Bcast
*** reported by process [1070923777,0]
*** on communicator MPI COMMUNICATOR 4 CREATE FROM 0
*** MPI_ERR_TYPE: invalid datatype
*** MPI_ERRORS_ARE_FATAL (processes in this communicator will now abort,
***    and potentially your MPI job)

Not sure if there is a mistake with the input files or if this is rather an installation problem of Intel oneAPI MKL and MPI. Here the the version of Open MPI I have installed (command line: $ ompi_info):

Package: Debian OpenMPI
Open MPI: 4.1.2

In the following the content of the structure.dat, input.dat and 4_h2o.in:

structure.dat:

12
ANGSTROM
O 1.6759513692647068 -0.6002853757772342 -1.6309012548029673 
H 1.0279586472537519 -0.9024715781296282 -2.26876660715738 
H 1.6152391970432078 -1.2290779155851816 -0.9104727361895202
O -1.5834963228239756 2.008349069216274 1.232856299827159
H -1.5747845721704647 1.0714599470204051 1.4334707422714463 
H -0.8617765374263848 2.128261335924955 0.6141278649996753 
O -1.3650917294404066 -0.565485109588074 2.1821923460132715
H -1.980150002454046 -0.467032993221947 -2.3069381847013353
H -1.350079743262641 -1.5072050335791736 2.00606312800933171
O 0.8016958070392737 2.1437494450214936 -0.1479543311321656
H 0.6999897857190609 1.6465013944331024 -0.9606557140636844
H 1.6273193428536925 1.830277846952749 0.22375067383692304

input.dat:

<method>
PIMD

<ensemble>
NVT

<bath_type>
MNHC

<temperature>
300.d0

<time_bath>
10.d0

<nnhc>
4

<dt>
0.25d0

<nstep>
10000

<corrections>
0  0

<nbead>
2

<ipotential>
QE

<np_beads>
2

<qe_input_file_name>
4_h2o.in

<qe_output>
1 10

<iboundary>
0

<iread_exit>
10

<iprint_std>
1

<iprint_bond>
10

<iprint_eavg>
10

<iprint_rest>
10

<iprint_trj>
10

<iprint_xyz>
1

<iprint_dip>
100

<iprint_rdf>
10

4_h2o.in:

&control
    calculation = 'scf'
    restart_mode='from_scratch'
    prefix = '4_h2o'
    outdir='./outdir'
    pseudo_dir = '../data/pseudopotentials'
/
&system
    ibrav=1, celldm(1) =9.86, nat=12, ntyp=2,
    ecutwfc=25, nosym=.true.
/
&electrons
    conv_thr=1e-8
/
&ions
/
ATOMIC_SPECIES
 H  1.00784  H.pbe-rrkjus.UPF
 O  15.999   O.pbe-rrkjus.UPF
ATOMIC_POSITIONS (angstrom)
O 1.6759513692647068 -0.6002853757772342 -1.6309012548029673 
H  1.0279586472537519 -0.9024715781296282 -2.26876660715738 
H  1.6152391970432078 -1.2290779155851816 -0.9104727361895202
O -1.5834963228239756 2.008349069216274 1.232856299827159
H  -1.5747845721704647 1.0714599470204051 1.4334707422714463 
H  -0.8617765374263848 2.128261335924955 0.6141278649996753 
O -1.3650917294404066 -0.565485109588074 2.1821923460132715
H  -1.980150002454046 -0.467032993221947 -2.3069381847013353
H  -1.350079743262641 -1.5072050335791736 2.00606312800933171
O 0.8016958070392737 2.1437494450214936 -0.1479543311321656
H  0.6999897857190609 1.6465013944331024 -0.9606557140636844
H  1.6273193428536925 1.830277846952749 0.22375067383692304 
K_POINTS (automatic)
  1 1 1 0 0 0

[MQS-mark]

I also built again PIMD and here the output. The cmake warning below could also be the reason for the above error (this was actually the problem and now I have gotten one step further. I re-built PIMD not in a Python environment to prevent libiomp5.so being set to the path of the MKL files for the Python environment):

cmake -DQE=on -DQEVERSION=6.3 -DAENET=on -DMKL=on ..
CMake Warning (dev) at /home/mark/miniconda3/envs/cebule/lib/python3.12/site-packages/cmake/data/share/cmake-3.30/Modules/ExternalProject/shared_internal_commands.cmake:1282 (message):
  The DOWNLOAD_EXTRACT_TIMESTAMP option was not given and policy CMP0135 is
  not set.  The policy's OLD behavior will be used.  When using a URL
  download, the timestamps of extracted files should preferably be that of
  the time of extraction, otherwise code that depends on the extracted
  contents might not be rebuilt if the URL changes.  The OLD behavior
  preserves the timestamps from the archive instead, but this is usually not
  what you want.  Update your project to the NEW behavior or specify the
  DOWNLOAD_EXTRACT_TIMESTAMP option with a value of true to avoid this
  robustness issue.
Call Stack (most recent call first):
  /home/mark/miniconda3/envs/cebule/lib/python3.12/site-packages/cmake/data/share/cmake-3.30/Modules/ExternalProject.cmake:3035 (_ep_add_download_command)
  lib/qe/CMakeLists.txt:35 (ExternalProject_Add)
This warning is for project developers.  Use -Wno-dev to suppress it.

-- URL=http://ann.atomistic.net/files/aenet-2.0.3.tar.bz2
-- aenet-2.0.3.tar.bz2 already downloaded.
-- MKL_ARCH: intel64
-- MKL_LINK: dynamic
-- MKL_INTERFACE_FULL: gf_lp64
-- MKL_THREADING: intel_thread
-- MKL_MPI: intelmpi
-- Found MKL: intelmpi
-- Found MKL: /opt/intel/oneapi/mkl/2024.2
-- Found MKL: /usr/lib/x86_64-linux-gnu/libmkl_gf_lp64.so
-- Found MKL: /usr/lib/x86_64-linux-gnu/libmkl_intel_thread.so
-- Found MKL: /usr/lib/x86_64-linux-gnu/libmkl_core.so
-- Found MKL: /home/mark/miniconda3/envs/cebule/lib/libiomp5.so
-- Configuring done (1.6s)
CMake Warning at CMakeLists.txt:47 (add_executable):
  Cannot generate a safe runtime search path for target pimd.mpi.x because
  files in some directories may conflict with libraries in implicit
  directories:

    runtime library [libz.so.1] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_gf_lp64.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_intel_thread.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_core.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib

  Some of these libraries may not be found correctly.

CMake Warning at CMakeLists.txt:48 (add_executable):
  Cannot generate a safe runtime search path for target pimd.x because files
  in some directories may conflict with libraries in implicit directories:

    runtime library [libmkl_gf_lp64.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_intel_thread.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_core.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib

  Some of these libraries may not be found correctly.

CMake Warning at CMakeLists.txt:49 (add_executable):
  Cannot generate a safe runtime search path for target polymers.x because
  files in some directories may conflict with libraries in implicit
  directories:

    runtime library [libz.so.1] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_gf_lp64.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_intel_thread.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib
    runtime library [libmkl_core.so] in /usr/lib/x86_64-linux-gnu may be hidden by files in:
      /home/mark/miniconda3/envs/cebule/lib

  Some of these libraries may not be found correctly.

[MQS-mark]

The installation/setup of PIMD+QE is now fixed (see previous comment).

The Si example also runs successfully.

For the 4 water molecules example I receive the following error:

Ensemble: NVT.

Temperature is controlled using massive Nose-Hoover chains thermostat.

System:    12 atoms,    4 beads,    300.0 K.

------------------------------------------------------------------------------
     type  symbol  number           mass            atoms
------------------------------------------------------------------------------
        1  O            8       15.99491                4
        2  H            1        1.00783                8
------------------------------------------------------------------------------

Atomic positions read from structure.dat. Molecular configuration initialized.

Free boundary condition - isolated system.

Thermostat positions and velocities initialized, given randomly.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine coset (1):
     nsym == 0
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     stopping ...
Abort(1) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0

In the following the input files:

input.dat

<input_style>
NEW

<method>
PIMD

<fdiff>
1.d-3

<nbead>
4

<ensemble>
NVT

<bath_type>
MNHC

<temperature>
300.d0

<nnhc>
4

<dt>
0.25d0

<nstep>
1000

<ipotential>
QE

<qe_input_file_name>
4_h2o.in

<qe_output>
1 10

<iboundary>
0

<iread_exit>
10

<iprint_std>
1

<iprint_bond>
10

<iprint_eavg>
10

<iprint_rest>
10

<iprint_trj>
10

<iformat_trj>
2

<iprint_xyz>
1

<iprint_dip>
100

<iprint_rdf>
10

<iprint_rgy>
10

<iprint_box>
10

4_h2o.in:

&control
    calculation = 'md',
    restart_mode='from_scratch',
    prefix = '4_h2o',
    outdir='./qe_outdir',
    dt=20,
    nstep=100,
    disk_io='high',
    iprint= 30,
    pseudo_dir = '../../pseudopotentials'
/
&system
    ibrav=1, celldm(1)=12.0, nat=12, ntyp=2,
    ecutwfc=25, ecutrho = 800.0, nosym=.true.
/
&electrons
    conv_thr=1.0e-8
/
&ions
/
ATOMIC_SPECIES
 H  1.00784  H.pbe-rrkjus.UPF
 O  15.999   O.pbe-rrkjus.UPF
ATOMIC_POSITIONS (angstrom)
O 1.6759513692647068 -0.6002853757772342 -1.6309012548029673 
H  1.0279586472537519 -0.9024715781296282 -2.26876660715738 
H  1.6152391970432078 -1.2290779155851816 -0.9104727361895202
O -1.5834963228239756 2.008349069216274 1.232856299827159
H  -1.5747845721704647 1.0714599470204051 1.4334707422714463 
H  -0.8617765374263848 2.128261335924955 0.6141278649996753 
O -1.3650917294404066 -0.565485109588074 2.1821923460132715
H  -1.980150002454046 -0.467032993221947 -2.3069381847013353
H  -1.350079743262641 -1.5072050335791736 2.00606312800933171
O 0.8016958070392737 2.1437494450214936 -0.1479543311321656
H  0.6999897857190609 1.6465013944331024 -0.9606557140636844
H  1.6273193428536925 1.830277846952749 0.22375067383692304 
K_POINTS (automatic)
  1 1 1 0 0 0

[MQS-mark]

@cometscome

The underlying problem of the above error message:

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine coset (1):
     nsym == 0
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

seems the lattice definition to be. I see the following notifications in the QE output log:

Message from routine volume:
     strange lattice parameter

    --------------------------------------------------------------------------------
NOTICE: Bravais lattice has wrong number (32) of symmetries - symmetries are disabled
--------------------------------------------------------------------------------
     Message from routine set_sym_bl:
     NOTICE: Symmetry group for Bravais lattice is not a group (2) - symmetries are disabled

And then there is no further output in the log file after this last notice. Seems to stop there. I just wonder what is wrong here because the example QE input file runs with QE without PIMD.

[cometscome]

Thank you for giving us this info. Now I asked Dr. Shiga. When I get something, I can tell you.

[MQS-mark]

@cometscome Thank you for all your help.

Here the log file of quantum espresso:

     Program PWSCF v.6.3MaX starts on 16Sep2024 at 16:59:39 

     This program is part of the open-source Quantum ESPRESSO suite
     for quantum simulation of materials; please cite
         "P. Giannozzi et al., J. Phys.:Condens. Matter 21 395502 (2009);
         "P. Giannozzi et al., J. Phys.:Condens. Matter 29 465901 (2017);
          URL http://www.quantum-espresso.org", 
     in publications or presentations arising from this work. More details at
     http://www.quantum-espresso.org/quote

     Parallel version (MPI), running on     1 processors

     MPI processes distributed on     1 nodes
     Reading input from 4_h2o.in

     Current dimensions of program PWSCF are:
     Max number of different atomic species (ntypx) = 10
     Max number of k-points (npk) =  40000
     Max angular momentum in pseudopotentials (lmaxx) =  3
               file H.pbe-rrkjus.UPF: wavefunction(s)  1S renormalized
               file O.pbe-rrkjus.UPF: wavefunction(s)  2S renormalized

     -- for PIMD -- 
     nstep in QE input            100
     nstep in PIMD input         1000
     nstep is forced to be the value of PIMD input
     In qe_init0 ::  nstep =         1000
     ------------- 

     -- flag for PIMD -- 
     lscf,lforce,lstres,lmd ::  T T F T
     ------------- 

     Message from routine volume:
     strange lattice parameter

     ***************************************
          PIMD step =          0
     ***************************************

     ***************************************
          alat [a.u.]         0.000000
          omega[a.u.^3]            NaN
     ***************************************

     ***************************************
          lattice vector [a.u.]
     ***************************************
       1         NaN         NaN         NaN
       2         NaN         NaN         NaN
       3         NaN         NaN         NaN
     ***************************************

     ***************************************
          lattice vector [-]
     ***************************************
       1         NaN         NaN         NaN
       2         NaN         NaN         NaN
       3         NaN         NaN         NaN
     ***************************************

     ***************************************
          coord_xyz [a.u.]
     ***************************************
       1    3.139919   -1.128887   -3.058755
       2    2.068280   -1.768979   -4.333447
       3    3.020137   -2.434044   -1.732297
       4   -2.943089    3.842499    2.373725
       5   -2.878365    1.925005    2.580987
       6   -1.634999    3.967683    0.979052
       7   -2.553457   -1.095532    4.097237
       8   -3.866135   -0.887641   -4.426357
       9   -2.325464   -2.615714    3.738713
      10    1.507485    4.140334   -0.261914
      11    1.334374    3.240556   -1.788474
      12    2.987467    3.406099    0.526840
     ***************************************

     ***************************************
          atomic position [-]
     ***************************************
       1    Infinity   -Infinity   -Infinity
       2    Infinity   -Infinity   -Infinity
       3    Infinity   -Infinity   -Infinity
       4   -Infinity    Infinity    Infinity
       5   -Infinity    Infinity    Infinity
       6   -Infinity    Infinity    Infinity
       7   -Infinity   -Infinity    Infinity
       8   -Infinity   -Infinity   -Infinity
       9   -Infinity   -Infinity    Infinity
      10    Infinity    Infinity   -Infinity
      11    Infinity    Infinity   -Infinity
      12    Infinity    Infinity    Infinity
     ***************************************

     gamma-point specific algorithms are used

     Subspace diagonalization in iterative solution of the eigenvalue problem:
     a serial algorithm will be used

--------------------------------------------------------------------------------
NOTICE: Bravais lattice has wrong number (32) of symmetries - symmetries are disabled
--------------------------------------------------------------------------------
     Message from routine set_sym_bl:
     NOTICE: Symmetry group for Bravais lattice is not a group (2) - symmetries are disabled

[cometscome]

I got the comment from Dr. Shiga. If you want to use QE, you have to set the boundary in input.dat, like

<iboundary>
  1
  9.232239523 -4.170626664  0.000000000
  0.000000000  7.223738307  0.000000000
  0.000000000  0.000000000  9.132346975

This example is for SiO2. please see examples/SiO2/qe_md.

[MQS-mark]

@cometscome Hi Yuki, This worked. The calculation runs through. I'm analyzing the results and will post all needed files for a 4 molecules water molecules system when the results are validated. Hope that this example is useful for your repository then.