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! module siesta_master
10!
11! Keeps a repository of coordinates, forces, and other properties, and handles
12! their transmission to/from a master program that uses siesta as a subroutine.
13! Since siesta was written as a main program, not as a subroutine, the flow of
14! data to/from the master program is not straighforward as subroutine arguments
15! and it was simpler to adress it with this module. Furthermore, there are two
16! communication modes available with the master program, through unix pipes and
17! through MPI, what justifies the following flow organization:
18! Using unix pipes (master+fsiesta_pipes compiled together):
19!   master <--> fsiesta_pipes <--pipe--> iopipes <--> siesta_master <--> siesta
20! Using MPI (master+siesta compiled together):
21!   master <--> fsiesta_mpi <--> siesta_master <--> siesta
22!
23!-----------------------------------------------------------------------------
24!
25! Public procedures provided by this module:
26!   coordsFromMaster     : Get atomic coordinates from master program
27!   forcesToMaster       : Send atomic forces to master program
28!   getForcesForMaster   : Get atomic forces to be sent to master program
29!   getPropertyForMaster : Get other properties to be sent to master program
30!   setPropertyForMaster : Set other properties to be sent to master program
31!   setCoordsFromMaster  : Set atomic positions sent by master program
32!   setMasterUnits       : Set physical units used by master program
33!
34! Public variables:
35!   siesta_server        : Is siesta a server?
36!   siesta_subroutine    : Is siesta a subroutine?
37!   input_file           : Name of fdf data file
38!
39!-----------------------------------------------------------------------------
40!
41! Interfaces of public procedures:
42!
43! subroutine coordsFromMaster( na, xa, cell )
44!   integer, intent(in) :: na        ! Number of atoms
45!   real(dp),intent(out):: xa(3,na)  ! Atomic coordinates
46!   real(dp),intent(out):: cell(3,3) ! Unit cell vectors
47! end subroutine coordsFromMaster
48!
49! subroutine forcesToMaster( na, Etot, fa, stress )
50!   integer, intent(in):: na          ! Number of atoms
51!   real(dp),intent(in):: Etot        ! Total energy
52!   real(dp),intent(in):: fa(3,na)    ! Atomic forces
53!   real(dp),intent(in):: stress(3,3) ! Stress tensor
54! end subroutine coordsFromMaster
55!
56! subroutine getForcesForMaster( na, fa, stress )
57!   integer, intent(in) :: na          ! Number of atoms
58!   real(dp),intent(out):: Etot        ! Total energy
59!   real(dp),intent(out):: fa(3,na)    ! Atomic forces
60!   real(dp),intent(out):: stress(3,3) ! Stress tensor
61! end subroutine getForcesForMaster
62!
63! subroutine setCoordsFromMaster( na, xa, cell )
64!   integer, intent(in):: na        ! Number of atoms
65!   real(dp),intent(in):: xa(3,na)  ! Atomic coordinates
66!   real(dp),intent(in):: cell(3,3) ! Unit cell vectors
67! end subroutine setCoordsFromMaster
68!
69! subroutine setMasterUnits( xunit, eunit )
70!   character(len=*),intent(in):: xunit  ! Physical unit of length
71!   character(len=*),intent(in):: eunit  ! Physical unit of energy
72! end subroutine setMasterUnits
73!
74! subroutine getPropertyForMaster( property, valueSize, value, units, error )
75!   character(len=*),intent(in) :: property   ! Property name
76!   integer,         intent(in) :: valueSize  ! Size of value array
77!   real(dp),        intent(out):: value(:)   ! Property value(s)
78!   character(len=*),intent(out):: units      ! Physical units
79!   character(len=*),intent(out):: error      ! Error name
80! end subroutine getPropertyForMaster
81!
82! subroutine setPropertyForMaster( property, valueSize, value, units )
83!   character(len=*),intent(in) :: property   ! Property name
84!   integer,         intent(in) :: valueSize  ! Size of value array
85!   real(dp),        intent(in) :: value(*)   ! Property value(s)
86!   character(len=*),intent(in) :: units      ! Physical units
87! end subroutine setPropertyForMaster
88!
89!-----------------------------------------------------------------------------
90
91MODULE siesta_master
92
93! Used module parameters and procedures
94  use precision, only: dp              ! Double precision real kind
95  use sys,       only: die             ! Termination routine
96  use fdf,       only: fdf_get         ! Reading fdf-options
97  use fdf,       only: fdf_convfac     ! Conversion of physical units
98
99  use iosockets, only: coordsFromSocket, forcesToSocket
100  use iopipes,   only: coordsFromPipe  ! Read coordinates from pipe
101  use iopipes,   only: forcesToPipe    ! Write forces to pipe
102
103  implicit none
104
105! Public procedures
106PUBLIC :: &
107  coordsFromMaster,     &! Get atomic coordinates from master program
108  forcesToMaster,       &! Send atomic forces to master program
109  getForcesForMaster,   &! Get atomic forces to be sent to master program
110  getPropertyForMaster, &! Get other properties to be sent to master program
111  setPropertyForMaster, &! Set other properties to be sent to master program
112  setCoordsFromMaster,  &! Set atomic positions sent by master program
113  setMasterUnits         ! Set physical units used by master program
114
115! Public variables
116  logical,public,save:: siesta_server     = .false. ! Is siesta a server?
117  logical,public,save:: siesta_subroutine = .false. ! Is siesta a subroutine?
118  character(len=132),public,save:: input_file = ' ' ! fdf data file
119
120PRIVATE ! Nothing is declared public beyond this point
121
122! Derived type to hold a physical property
123  type propType
124    private
125    character(len=32):: name = ' '
126    character(len=32):: units = ' '
127    integer          :: size = 0
128    real(dp),pointer :: value(:)
129  end type propType
130
131! Physical units used by siesta
132  character(len=*), parameter :: siesta_xunit = 'Bohr'
133  character(len=*), parameter :: siesta_eunit = 'Ry'
134  character(len=*), parameter :: siesta_funit = 'Ry/Bohr'
135  character(len=*), parameter :: siesta_sunit = 'Ry/Bohr**3'
136
137! Global module variables and default values
138  integer,     parameter:: maxProps = 100      ! Max number of properties
139  integer,          save:: nProps   = 0        ! Number of stored properties
140  integer,          save:: nAtoms   = 0        ! Number of atoms
141  real(dp),         save:: ucell(3,3) = 0      ! Unit cell vectors
142  real(dp),         save:: stressT(3,3) = 0    ! Stress tensor
143  real(dp),         save:: energy = 0          ! Total energy
144  real(dp),pointer, save:: coords(:,:)         ! Atomic coordinates
145  real(dp),pointer, save:: forces(:,:)         ! Atomic forces
146  type(propType),   save:: prop(maxProps)      ! Stored physical properties
147  character(len=32),save:: master_xunit = 'Ang'       ! Master's unit of length
148  character(len=32),save:: master_eunit = 'eV'        ! Master's unit of energy
149  character(len=32),save:: master_funit = 'eV/Ang'    ! Master's unit of force
150  character(len=32),save:: master_sunit = 'eV/Ang**3' ! Master's unit of stress
151
152CONTAINS
153
154!-----------------------------------------------------------------------------
155
156subroutine coordsFromMaster( na, xa, cell )
157
158  implicit none
159  integer, intent(in) :: na        ! Number of atoms
160  real(dp),intent(out):: xa(3,na)  ! Atomic coordinates
161  real(dp),intent(out):: cell(3,3) ! Unit cell vectors
162  character*32 :: iface            ! Interface mode
163
164! In the pipes version, this is where we first know that we are a server
165  siesta_server = .true.
166
167  if (siesta_subroutine) then
168    if (na==nAtoms) then
169      xa = coords
170      cell = ucell
171    else
172      call die('coordsFromMaster: ERROR: number-of-atoms mismatch')
173    endif
174  else
175    iface = fdf_get( "Master.interface",  "pipe")
176    if ( iface == "pipe") call coordsFromPipe( na, xa, cell )
177    if ( iface == "socket") call coordsFromSocket (na, xa, cell )
178  end if ! (siesta_subroutine)
179
180end subroutine coordsFromMaster
181
182!-----------------------------------------------------------------------------
183
184subroutine forcesToMaster( na, Etot, fa, stress )
185
186  implicit none
187  integer, intent(in):: na          ! Number of atoms
188  real(dp),intent(in):: Etot        ! Total energy
189  real(dp),intent(in):: fa(3,na)    ! Atomic forces
190  real(dp),intent(in):: stress(3,3) ! Stress tensor
191  character*32 :: iface            ! Interface mode
192
193  if (siesta_subroutine) then
194    if (na==nAtoms) then
195      energy = Etot
196      forces = fa
197      stressT = stress
198    else ! (na/=nAtoms)
199      call die('coordsFromMaster: ERROR: number-of-atoms mismatch')
200    endif ! (na==nAtoms)
201  else
202    iface = fdf_get( "Master.interface",  "pipe")
203    if ( iface == "pipe") call forcesToPipe( na, Etot, fa, stress )
204    if ( iface == "socket") call forcesToSocket( na, Etot, fa, stress )
205  end if ! (siesta_subroutine)
206
207end subroutine forcesToMaster
208
209!-----------------------------------------------------------------------------
210
211subroutine getForcesForMaster( na, Etot, fa, stress )
212
213  implicit none
214  integer, intent(in) :: na          ! Number of atoms
215  real(dp),intent(out):: Etot        ! Total energy
216  real(dp),intent(out):: fa(3,na)    ! Atomic forces
217  real(dp),intent(out):: stress(3,3) ! Stress tensor
218
219! Check that number of atoms is right, then copy data from repository
220! and convert them to master's physical units
221  if (na==nAtoms) then
222    Etot   = energy  * fdf_convfac( siesta_eunit, master_eunit )
223    fa     = forces  * fdf_convfac( siesta_funit, master_funit )
224    stress = stressT * fdf_convfac( siesta_sunit, master_sunit )
225  else
226    call die('getForcesForMaster: ERROR: number-of-atoms mismatch')
227  end if
228
229end subroutine getForcesForMaster
230
231!-----------------------------------------------------------------------------
232
233subroutine setCoordsFromMaster( na, xa, cell )
234
235  implicit none
236  integer, intent(in):: na        ! Number of atoms
237  real(dp),intent(in):: xa(3,na)  ! Atomic coordinates
238  real(dp),intent(in):: cell(3,3) ! Unit cell vectors
239
240! Allocate arrays for repository the first time
241  if (nAtoms==0) then
242    nAtoms = na
243    allocate( coords(3,na), forces(3,na) )
244  end if
245
246! Check that number of atoms is right, then copy data to repository
247! after converting them to siesta physical units
248  if (na==nAtoms) then
249    coords = xa   * fdf_convfac( master_xunit, siesta_xunit )
250    ucell  = cell * fdf_convfac( master_xunit, siesta_xunit )
251  else
252    call die('setCoordsFromMaster: ERROR: number-of-atoms mismatch')
253  end if
254
255end subroutine setCoordsFromMaster
256
257!-----------------------------------------------------------------------------
258
259subroutine setMasterUnits( xunit, eunit )
260
261  implicit none
262  character(len=*),intent(in):: xunit  ! Physical unit of length
263  character(len=*),intent(in):: eunit  ! Physical unit of energy
264
265  master_xunit = xunit
266  master_eunit = eunit
267  master_funit = trim(eunit)//'/'//trim(xunit)         ! Unit of force
268  master_sunit = trim(eunit)//'/'//trim(xunit)//'**3'  ! Unit of stress
269
270end subroutine setMasterUnits
271
272!-----------------------------------------------------------------------------
273
274subroutine getPropertyForMaster( property, valueSize, value, units, error )
275
276  implicit none
277  character(len=*),intent(in) :: property   ! Property name
278  integer,         intent(in) :: valueSize  ! Size of value array
279  real(dp),        intent(out):: value(:)   ! Property value(s)
280  character(len=*),intent(out):: units      ! Physical units
281  character(len=*),intent(out):: error      ! Error name
282
283  logical:: found
284  integer:: iProp
285
286  found = .false.
287  do iProp = 1,nProps
288    if (prop(iProp)%name == property) then
289      if (prop(iProp)%size == valueSize) then
290        value = prop(iProp)%value
291        units = prop(iProp)%units
292        error = ' '
293      else
294        error = 'wrong_size'
295      end if
296      found = .true.
297      exit ! iProp loop
298    end if
299  end do ! iProp
300  if (.not.found) error = 'unknown_property'
301
302end subroutine getPropertyForMaster
303
304!-----------------------------------------------------------------------------
305
306subroutine setPropertyForMaster( property, valueSize, value, units )
307
308  implicit none
309  character(len=*),intent(in) :: property   ! Property name
310  integer,         intent(in) :: valueSize  ! Size of value array
311  real(dp),        intent(in) :: value(*)   ! Property value(s)
312  character(len=*),intent(in) :: units      ! Physical units
313
314  logical:: found
315  integer:: iProp
316
317! Check if property is already stored
318  found = .false.
319  do iProp = 1,nProps
320    if (prop(iProp)%name == property) then
321      found = .true.
322      exit ! iProp loop
323    end if
324  end do ! iProp
325
326! Increase the number of stored properties, if necessary
327  if (.not.found) then
328    nProps = nProps + 1
329    iProp = nProps
330  end if
331
332! (Re)allocate the array to store value(s)
333  if (prop(iProp)%size /= valueSize) then
334    if (prop(iProp)%size /= 0) deallocate( prop(iProp)%value )
335    allocate( prop(iProp)%value(valueSize) )
336  end if
337
338! Store data
339  prop(iProp)%name  = property
340  prop(iProp)%size  = valueSize
341  prop(iProp)%value = value(1:valueSize)
342  prop(iProp)%units = units
343
344end subroutine setPropertyForMaster
345
346END MODULE siesta_master
347
348