1! ---
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt .
6! See Docs/Contributors.txt for a list of contributors.
7! ---
8
9!*****************************************************************************
10! module iosockets
11!*****************************************************************************
12! Handles communications between siesta-as-a-subroutine and a driver
13! program, transfering coordinates and forces trough POSIX sockets.
14! The routines that handle the other side of the communication are
15! in file fsiesta_sockets.f90 and/or in the i-PI code.
16! J.M.Soler, A.Garcia, and M.Ceriotti, 2014
17!*****************************************************************************
18!
19!   Modules used:
20! use precision,    only: dp
21! use parallel,     only: IOnode
22! use fdf
23! use f90sockets,   only: open_socket, writebuffer, readbuffer
24! use sys,          only: die, bye
25! use m_mpi_utils,  only: broadcast
26! use cellSubs,     only: volcel
27! use mpi_siesta
28!
29!   Public derived types:
30! none
31!
32!   Public variables and arrays:
33! none
34!
35!   Public procedures:
36! coordsFromSocket, &! receive coords and cell vectors from master program
37! forcesToSocket     ! send energy, forces, and stress to master program
38!
39!*****************************************************************************
40! subroutine coordsFromSocket( na, xa, cell )
41!   integer,  intent(in) :: na            ! Number of atoms
42!   real(dp), intent(out):: xa(3,na)      ! Atomic coordinates (bohr)
43!   real(dp), intent(out):: cell(3,3)     ! Lattice vectors (bohr)
44!*****************************************************************************
45! subroutine forcesToSocket( na, energy, forces, stress )
46!   integer,  intent(in) :: na            ! Number of atoms
47!   real(dp), intent(in) :: energy        ! Total energy (Ry)
48!   real(dp), intent(in) :: forces(3,na)  ! Atomic forces (Ry/bohr)
49!   real(dp), intent(in) :: stress(3,3)   ! Stress tensor (Ry/bohr^3)
50!*****************************************************************************
51
52module iosockets
53
54! Used modules
55  use precision,    only: dp
56  use parallel,     only: IOnode
57  use fdf
58  use f90sockets,   only: open_socket, writebuffer, readbuffer, close_socket
59  use sys,          only: die, bye
60  use m_mpi_utils,  only: broadcast
61  use cellSubs,     only: volcel
62#ifdef MPI
63  use mpi_siesta
64#endif
65
66  implicit none
67
68! Public procedures
69PUBLIC :: &
70  coordsFromSocket, &! receive coords and cell vectors from master program
71  forcesToSocket     ! send energy, forces, and stress to master program
72
73PRIVATE  ! Nothing is declared public beyond this point
74
75! Private module parameters
76! WARNING: IPI_MSGLEN and MSGLEN must be equal in i-PI and fsiesta, respectively
77  integer,                 parameter ::   IPI_MSGLEN = 12
78  integer,                 parameter ::       MSGLEN = 80
79  integer,                 parameter ::     unit_len = 32
80  character(len=unit_len), parameter :: siesta_xunit = 'bohr'
81  character(len=unit_len), parameter :: siesta_eunit = 'ry'
82  character(len=unit_len), parameter ::    ipi_xunit = 'bohr'
83  character(len=unit_len), parameter ::    ipi_eunit = 'hartree'
84
85! Private module variables
86  integer, save :: socket
87  real(dp),save :: cellv
88  character(len=32),      save:: master='unknown', stype='unknown'
89  character(len=unit_len),save:: master_eunit='unknown', &
90                                 master_xunit='unknown'
91#ifdef MPI
92  integer :: MPIerror
93#endif
94
95CONTAINS
96
97!*****************************************************************************
98
99subroutine coordsFromSocket( na, xa, cell )
100  ! Reads coordinates from socket
101  implicit none
102  integer,  intent(in)  :: na         ! Number of atoms
103  real(dp), intent(out) :: xa(3,na)   ! Atomic coordinates (bohr)
104  real(dp), intent(out) :: cell(3,3)  ! Lattice vectors (bohr)
105
106! Local parameters
107  character(len=*),parameter:: myName='coordsFromSocket '
108
109! Local variables and arrays
110  logical, save            :: firstTime = .true.
111  CHARACTER(len=MSGLEN)    :: header, message
112  CHARACTER(len=1024)      :: host
113  INTEGER                  :: inet, port, n
114  REAL(dp)                 :: aux(9), c(9)
115  REAL(dp),ALLOCATABLE     :: x(:)
116
117! Open the socket, shared by the receive and send sides of communication
118  if (firstTime) then
119    master = fdf_get( "Master.code",  "fsiesta")
120    host = fdf_get(   "Master.address",  "localhost")
121    port = fdf_get(   "Master.port",  10001)
122    stype = fdf_get(  "Master.socketType",  "inet")
123
124    if (leqi(stype,'unix')) then
125      inet = 0
126    elseif (leqi(stype,'inet')) then
127      inet = 1
128    else
129      call die(myName//'ERROR: unknown socket type '//trim(stype))
130    endif
131
132    if (IOnode) then
133      write(*,'(/,a,i4,i8,2x,a)') &
134          myName//'opening socket: inet,port,host=',inet,port,trim(host)
135      call open_socket(socket, inet, port, host)
136    endif
137
138    firstTime = .false.
139  end if ! (firstTime .and. IOnode)
140
141! Read header from socket
142  header = ''
143  if (IOnode) then
144    do
145      if (leqi(master,'i-pi')) then
146        call readbuffer(socket, header, IPI_MSGLEN)
147
148        ! Immediately stop if requested
149        if (trim(header)=='EXIT') then ! we are done!
150          call close_socket(socket)
151          call bye(myName//'STOP requested by driver')
152        end if
153
154        if (trim(header)/='STATUS') exit ! do loop
155        message = 'READY'
156        call writebuffer(socket, message, IPI_MSGLEN)
157      elseif (leqi(master,'fsiesta')) then
158        call readbuffer(socket, header)
159
160        ! Immediately stop if requested
161        if (trim(header)=='quit') then
162          call writebuffer(socket,'quitting')
163          call close_socket(socket)
164          call bye(myName//'STOP requested by driver')
165        elseif (trim(header)=='wait') then
166          call writebuffer(socket,'ready')
167          cycle ! do loop
168        else
169          exit ! do loop
170        endif
171      else
172        call die(myName//'ERROR: unknown master')
173      endif ! leqi(master)
174    enddo
175  endif ! IOnode
176#ifdef MPI
177  call MPI_Bcast(header, MSGLEN, MPI_Character, 0, MPI_Comm_World, MPIerror)
178#endif
179
180! Read cell vectors from socket, as a single buffer vector
181  if (IOnode) then
182    if (leqi(master,"i-pi") .and. trim(header)=="POSDATA") then
183      call readbuffer(socket, c, 9)
184      ! Please see below for the cell array
185      ! If this is to be used, it should *ALSO* be transposed!
186      call readbuffer(socket, aux, 9)    ! not used in siesta
187      master_xunit = ipi_xunit
188      master_eunit = ipi_eunit
189    elseif (leqi(master,"fsiesta") .and. trim(header)=="begin_coords") then
190      call readbuffer(socket, master_xunit)
191      call readbuffer(socket, master_eunit)
192      call readbuffer(socket, c, 9)
193    else
194      call die(myName//'ERROR: unexpected header: '//trim(header))
195    end if
196  endif
197
198! Broadcast and copy cell from buffer
199#ifdef MPI
200  call MPI_Bcast(master_xunit, unit_len, MPI_Character,0,  &
201                 MPI_Comm_World, MPIerror)
202  call MPI_Bcast(master_eunit, unit_len, MPI_Character,0,  &
203                 MPI_Comm_World, MPIerror)
204  call MPI_Bcast(c,9, MPI_Double_Precision,0, MPI_Comm_World, MPIerror)
205#endif
206  ! I-Pi assumes row-major order of *only* the cell parameters
207  ! So we need to transpose the cell
208  ! This is a very bad practice since the same is happening in
209  ! ASE I-pi implementation.
210  ! The cell is transposed when transfering, but not the coordinates.
211  ! Essentially one could leave out *any* transposes in all implementations
212  ! and it would work! However, this is now a legacy implementation
213  ! that should not be changed, neither in I-pi, nor ASE.
214  cell = TRANSPOSE( RESHAPE( c, (/3,3/) ) )
215
216! Read and check number of atoms
217  if (ionode) then
218    call readbuffer(socket, n)
219    if (n/=na) call die(myName//'ERROR: unexpected number of atoms')
220  endif
221
222! Read and broadcast atomic coordinates
223  allocate(x(3*na))
224  if (IOnode) call readbuffer(socket, x, 3*na)
225#ifdef MPI
226  call MPI_Bcast(x,3*na, MPI_Double_Precision,0, MPI_Comm_World, MPIerror)
227#endif
228  xa = RESHAPE( x, (/3,na/) )
229  deallocate(x)
230
231! Read trailing message
232  if (IOnode .and. leqi(master,'fsiesta')) then
233    call readbuffer(socket, message)
234    if (message/='end_coords') then
235      call die(myName//'ERROR: unexpected trailer:'//trim(message))
236    end if
237  endif
238
239  if ( IONode ) then
240    ! Print coordinates and cell vectors received
241    write(*,'(/,4a,/,(3f12.6))') myName,'cell (',trim(master_xunit),') =', cell
242    write(*,'(4a,/,(3f12.6))') myName,'coords (',trim(master_xunit),') =', xa
243  end if
244
245! Convert physical units
246  cell = cell * fdf_convfac( master_xunit, siesta_xunit )
247  xa   = xa   * fdf_convfac( master_xunit, siesta_xunit )
248
249! Compute and save the cell volume, to be used by forcesToSocket
250  cellv = volcel(cell)
251
252end subroutine coordsFromSocket
253
254!*****************************************************************************
255
256subroutine forcesToSocket( na, energy, forces, stress )
257! Writes stress and forces to socket
258  implicit none
259  integer,  intent(in) :: na            ! Number of atoms
260  real(dp), intent(in) :: energy        ! Total energy (Ry)
261  real(dp), intent(in) :: forces(3,na)  ! Atomic forces (Ry/bohr)
262  real(dp), intent(in) :: stress(3,3)   ! Stress tensor (Ry/bohr^3)
263
264! Local parameters
265  character(len=*),parameter:: myName='forcesToSocket '
266
267! Local variables and arrays
268  character(len=MSGLEN):: header, message
269  real(dp)             :: e, s(9), vir(9)
270  real(dp),allocatable :: f(:)
271
272! Copy input to local variables
273  allocate(f(3*na))
274  e = energy
275  f(:) = reshape( forces, (/3*na/) )
276  s(:) = reshape( stress, (/9/) )
277
278! Convert physical units
279  e = e * fdf_convfac( siesta_eunit, master_eunit )
280  f(:) = f(:) * fdf_convfac( siesta_eunit, master_eunit ) &
281        / fdf_convfac( siesta_xunit, master_xunit )
282  s(:) = s(:) * fdf_convfac( siesta_eunit, master_eunit ) &
283        / fdf_convfac( siesta_xunit, master_xunit )**3
284
285! Find virial tensor for i-pi, in master's units
286  vir(:) = -s * cellv * fdf_convfac( siesta_xunit, master_xunit )**3
287
288! Write forces to socket
289  if (IOnode) then
290    if (leqi(master,'i-pi')) then
291      do
292        call readbuffer(socket, header, IPI_MSGLEN)
293        if (trim(header)=='STATUS') then        ! inform i-pi of my status
294          message = 'HAVEDATA'
295          call writebuffer(socket,message,IPI_MSGLEN)
296          cycle ! do loop
297        elseif (trim(header)=='GETFORCE') then  ! proceed to send forces
298          exit ! do loop
299        else
300          call die(myName//'ERROR: unexpected header from i-pi')
301        endif
302      enddo
303      message = 'FORCEREADY'
304      call writebuffer(socket,message,IPI_MSGLEN)
305      call writebuffer(socket,e)
306      call writebuffer(socket,na)
307      call writebuffer(socket,f,3*na)
308      call writebuffer(socket,vir,9)
309      ! i-pi can also receive an arbitrary string, that will be printed
310      ! out to the "extra" trajectory file. This is useful if you want
311      ! to return additional information, e.g. atomic charges, wannier
312      ! centres, etc. one must return the number of characters, then
313      ! the string. here we just send back zero characters.
314      call writebuffer(socket,0)
315    elseif (leqi(master,'fsiesta')) then
316!      do
317!        call readbuffer(socket, header)
318!        if (trim(header)/='wait') exit
319!      enddo
320      call writebuffer(socket,'begin_forces')
321      call writebuffer(socket,e)
322      call writebuffer(socket,s,9)
323      call writebuffer(socket,na)
324      call writebuffer(socket,f,3*na)
325      call writebuffer(socket,'end_forces')
326    else
327      call die(myName//'ERROR: unknown master')
328    endif ! leqi(master)
329  end if ! IOnode
330
331  if ( IONode ) then
332    ! Print energy, forces, and stress tensor sent to master
333    write(*,'(/,a,f12.6)') myName// 'energy ('//trim(master_eunit)//') =', e
334    write(*,'(a,/,(3f12.6))') myName//'stress ('//trim(master_eunit)//'/'//trim(master_xunit)//'^3) =', s
335    write(*,'(a,/,(3f12.6))') myName//'forces ('//trim(master_eunit)//'/'//trim(master_xunit)//') =', f
336  end if
337  deallocate(f)
338
339end subroutine forcesToSocket
340
341end module iosockets
342