Hi Aniruddha M Dive,
Below is the code that I changed a few lines in cp2k_shell.F and ase-cp2k calculator. ( You will see the lines I commented out)
Please note that the code is only tested wit cp2k-v6.1 (I found the most recent cp2k code revised cp2k_shell.F significantly).
The idea is that the python code will write a file "
CP2K_POSITIONS" with coordinates instead of writing them into PIPE,
and cp2k_shell will read the positions from the file instead of the PIPE.
Best,
Geng
IF (para_env%mepos==para_env%source) THEN
!READ (*,*,iostat=iostat) n_atom2
!IF (iostat/=0) CPABORT('setpos read n_atom')
!IF (n_atom2/=SIZE(pos)) THEN
! CALL my_assert(.FALSE.,'setpos invalid number of atoms',failure)
! DO i=1,n_atom
! READ(*,'(a)',iostat=iostat) cmdStr
! CALL compress(cmdStr,full=.TRUE.)
! CALL uppercase(cmdStr)
! IF (cmdStr=='*END') EXIT
! END DO
! GOTO 10
!END IF
!READ (*,*,iostat=iostat) pos
!IF (iostat/=0) CPABORT('setpos read coord')
inquire(unit=201,opened=unitalive)
if (unitalive) CPABORT('UNIT 201 is being used')
inquire(file="CP2K_POSITIONS",exist=filepresence)
if (.not. filepresence) CPABORT('FILE CP2K_POSITIONS NOT EXIST')
open(201,action='READ',file="CP2K_POSITIONS",iostat=iostat,status='OLD',form='FORMATTED',ACCESS='SEQUENTIAL')
if (iostat/=0) CPABORT('read CP2K_POSITIONS')
READ(201,*,iostat=iostat) n_atom2
IF (iostat/=0) CPABORT('setpos read n_atom2')
IF (n_atom2/=SIZE(pos)) THEN
CALL my_assert(.FALSE.,'setpos invalid number of atoms',failure)
DO i=1,n_atom
READ(201,'(a)',iostat=iostat) cmdStr
CALL compress(cmdStr,full=.TRUE.)
CALL uppercase(cmdStr)
IF (cmdStr=='*END') EXIT
END DO
GOTO 10
END IF
READ (201,*,iostat=iostat) pos
CLOSE(201)
IF (iostat/=0) CPABORT('setpos read coord')
pos(:) = pos(:)/pos_fact
READ(*,'(a)',iostat=iostat) cmdStr
CALL compress(cmdStr,full=.TRUE.)
CALL uppercase(cmdStr)
CALL my_assert(cmdStr=='*END',' missing *END',failure)
END IF
---changes in ase cp2k calcualtor -----
if 'positions' in system_changes:
with open("CP2K_POSITIONS","w") as fp:
fp.write("%d\n" % (3*n_atoms))
for pos in self.atoms.get_positions():
fp.write('%.18e %.18e %.18e\n' % tuple(pos))
self._shell.send('SET_POS %d' % self._force_env_id)
#self._shell.send('%d' % (3 * n_atoms))
#for pos in self.atoms.get_positions():
# self._shell.send('%.18e %.18e %.18e' % tuple(pos))
self._shell.send('*END')
max_change = float(self._shell.recv())
assert max_change >= 0 # sanity check
self._shell.expect('* READY')
if os.path.isfile("CP2K_POSITIONS"):
os.remove("CP2K_POSITIONS")