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! ---
8module jobList
9
10  implicit none
11
12PUBLIC:: &
13  countJobs,  &! count number of jobs and lists
14  runJobs,    &! runs a job list read from a file unit
15  getResults   ! collect requested results into summary files
16
17PRIVATE    ! nothing is declared public beyond this point
18
19  ! Internal parameters
20  character(len=*),parameter:: dataSuffix = '.fdf'
21  character(len=*),parameter:: endSuffix = '.EIG'
22  character(len=*),parameter:: file0Exit = '0_NORMAL_EXIT'
23  character(len=*),parameter:: defaultFiles = &
24                                       '*.fdf *.vps *.psf *.ion queue.sh'
25  character(len=*),parameter:: defaultQueue = &
26                                       './siesta < $jobName.fdf > $jobName.out'
27  character(len=*),parameter:: defaultRequest(1) = (/'energy'/)
28  integer,parameter:: maxLines = 1000   ! Max lines in job-list file
29  integer,parameter:: maxWords = 100    ! Max words in one line
30  integer,parameter:: maxJobs = 1000    ! Max jobs in job list
31  integer,parameter:: ll      = 500     ! Max characters per line
32  integer,parameter:: wl      = 500     ! Max characters per word
33  integer,parameter:: unitIn = 42       ! I/O unit for input files
34  integer,parameter:: unitOut = 43      ! I/O unit for output files
35  integer,parameter:: dp = kind(1.d0)
36
37  ! Internal variables and arrays
38  integer,save:: jobCores, totCores, totJobs, totLists
39
40  ! Derived type to hold calculation results
41  type resultsType
42    private
43    real(dp)        :: cell(3,3)   =0._dp
44    real(dp)        :: stress(3,3) =0._dp
45    real(dp)        :: volume      =0._dp
46    real(dp)        :: energy      =0._dp
47    real(dp)        :: pressure    =0._dp
48    real(dp)        :: virial      =0._dp
49    real(dp),pointer:: coord(:,:)  =>null()
50    real(dp),pointer:: force(:,:)  =>null()
51    integer, pointer:: za(:)       =>null()   ! atomic numbers
52  end type resultsType
53
54CONTAINS
55
56!------------------------------------------------------------------------------
57
58subroutine countJobs( unit, nLists, nJobs, nCores )
59
60! Counts number of jobs and lists in a joblist file, specified by a I/O unit
61
62  implicit none
63  integer,         intent(in) :: unit   ! IO unit of datafile
64  integer,optional,intent(out):: nLists ! total number of job lists
65  integer,optional,intent(out):: nJobs  ! total number of jobs
66  integer,optional,intent(out):: nCores ! total number of required cores
67
68  integer :: nc
69  character(len=wl):: myDir
70
71  ! Get current directory and add trailing '/' if necessary
72  call getcwd(myDir)
73  nc = len(trim(myDir))
74  if (myDir(nc:nc) /= '/') myDir = myDir(1:nc) // '/'
75
76  ! Scan the 'root' list (the whole file), specified by blank name
77  call scanList(unit,myDir,defaultQueue,defaultFiles,defaultRequest,' ','count')
78
79  ! Copy number of jobs and lists to output variables
80  if (present(nLists)) nLists = totLists
81  if (present(nJobs))  nJobs  = totJobs
82  if (present(nCores)) nCores = totCores
83
84end subroutine countJobs
85
86!------------------------------------------------------------------------------
87
88subroutine runJobs( unit )
89
90! Runs or queues a number of jobs in a joblist file, specified by a I/O unit
91
92  implicit none
93  integer,intent(in) :: unit   ! IO unit of datafile
94
95  integer :: nc
96  character(len=wl):: myDir
97
98  ! Get current directory and add trailing '/' if necessary
99  call getcwd(myDir)
100  nc = len(trim(myDir))
101  if (myDir(nc:nc) /= '/') myDir = myDir(1:nc) // '/'
102
103  ! Scan the 'root' list (the whole file), specified by blank name
104  call scanList(unit,myDir,defaultQueue,defaultFiles,defaultRequest,' ','run')
105
106end subroutine runJobs
107
108!------------------------------------------------------------------------------
109
110subroutine getResults( unit )
111
112! Collects requested results of jobs in a joblist file (specified by a
113! I/O unit), into summary output files, within each sub-list directory
114
115  implicit none
116  integer,intent(in) :: unit   ! IO unit of datafile
117
118  integer :: nc
119  character(len=wl):: myDir
120
121  ! Get current directory and add trailing '/' if necessary
122  call getcwd(myDir)
123  nc = len(trim(myDir))
124  if (myDir(nc:nc) /= '/') myDir = myDir(1:nc) // '/'
125
126  ! Scan the 'root' list (the whole file), specified by blank name
127  call scanList(unit,myDir,defaultQueue,defaultFiles,defaultRequest,' ','get')
128
129end subroutine getResults
130
131!------------------------------------------------------------------------------
132
133recursive subroutine scanList( unit, dir, queue, files, request, &
134                               listName, task )
135
136  implicit none
137  integer,         intent(in) :: unit       ! I/O unit of joblist file
138  character(len=*),intent(in) :: dir        ! parent directory
139  character(len=*),intent(in) :: queue      ! queuing statement
140  character(len=*),intent(in) :: files      ! required data files
141  character(len=*),intent(in) :: request(:) ! requested results
142  character(len=*),intent(in) :: listName   ! name of job list
143  character(len=*),intent(in) :: task       ! ('count'|'run'|'get')
144
145  character(len=1),parameter:: separator(1) = (/' '/)
146  integer :: iCase, iLine, iostat, nCases, nCores, nResults, nWords, nJobs
147  character(len=wl):: caseDir(maxJobs), caseName(maxJobs), caseType(maxJobs), &
148                       fileIn, fileOut, line, myDir, myFiles, myQueue, &
149                       myRequest(maxWords), newList, words(maxWords)
150  real(dp):: results(maxWords)
151  logical :: finished
152
153  ! Set the list directory and copy data files to it
154  if (listName==' ') then  ! this is the 'root' list (the whole file)
155    myDir = dir
156  else                     ! this is a genuine list of jobs
157    myDir = trim(dir) // trim(listName) // '/'
158    ! Create a new directory for this list and copy files to it
159    if (task=='run') then
160      call system('mkdir -p ' // trim(myDir))
161      call system('cp -f ' // trim(files) // ' ' // trim(myDir) // &
162                  ' 2> /dev/null' )
163      call chdir(trim(myDir))
164    endif
165  endif ! (listName==' ')
166
167  ! Global initializations
168  if (listName==' ') then   ! this is the 'root' list (the whole file)
169    totCores = 0            ! total number of cores
170    totJobs = 0             ! total number of jobs
171    totLists = 0            ! total number of job lists
172  endif
173
174  ! Set my own copies of queue, files, and request
175  nResults = count(request/=' ')
176  myRequest = ' '
177  myRequest(1:nResults) = request(1:nResults)
178  myQueue = queue
179  myFiles = files
180  call getCores(myQueue,jobCores)   ! get number of cores per job
181
182  ! Loop on lines of datafile
183  finished = .false.       ! has the list terminated normally?
184  nJobs = 0                ! number of jobs in list
185  nCores = 0               ! number of cores required by jobs in list
186  nCases = 0               ! number of 'cases' (jobs or sublists) in list
187  do iLine = 1,maxLines
188
189    ! Read one line (possibly continued) of input file, ignoring comments
190    call readLine( unit, line, iostat )
191
192    ! End of file check
193    if (iostat<0) then                      ! end of file
194      if (listName==' ') then               ! this is the 'root' list
195        finished = .true.
196        exit ! iLine loop
197      else
198        print*,'runList ERROR: unterminated %list = '//trim(listName)
199        stop
200      endif ! (listName==' ')
201    endif ! (iostat<0)
202
203    ! Parse line
204    call parser(separator,line,words,nWords)
205
206    ! Act depending on line content
207    if (nWords==0) then                     ! blank or comment line
208      cycle
209    elseif (trim(words(1))=='%queue') then  ! queue specification
210      line = adjustl(line)
211      myQueue = adjustl(line(7:))           ! remove '%queue' from line
212      call getCores(myQueue,jobCores)
213    elseif (trim(words(1))=='%files') then  ! data files specification
214      line = adjustl(line)
215      myFiles = adjustl(line(7:))           ! remove '%files' from line
216    elseif (trim(words(1))=='%result') then ! result-request especification
217      myRequest = ' '
218      myRequest(1:nWords-1) = words(2:nWords) ! exclude 1st word (it is %result)
219    elseif (words(1)=='%list') then         ! new sub-list
220      nCases = nCases+1
221      newList = words(2)
222      caseType(nCases) = 'list'
223      caseName(nCases) = newList
224      caseDir(nCases) = trim(myDir) // trim(newList) // '/'
225      call scanList(unit,myDir,myQueue,myFiles,myRequest,newList,task)
226    elseif (words(1)=='%endlist') then      ! end of my list
227      if (words(2)==listName) then
228        finished = .true.
229        exit ! iLine loop
230      else
231        print*,'runList ERROR: inconsistent %list = '//trim(listName)// &
232               ', %endList = '//trim(words(2))
233        stop
234      endif
235    else                                    ! new job
236      nCases = nCases+1
237      nJobs = nJobs+1
238      nCores = nCores + jobCores
239      caseType(nCases) = 'job'
240      call nameJob( line, caseName(nCases) )
241      caseDir(nCases) = trim(myDir) // trim(caseName(nCases)) // '/'
242      if (task=='run') call runOneJob(myDir,myQueue,myFiles,line)
243    endif
244
245  end do ! iLine
246  if (.not.finished) stop 'runList ERROR: parameter maxLines too small'
247
248  ! Write file of requested results
249  if (task=='get') then
250    if (listName==' ') then
251      fileOut = trim(myDir)//'jobList.results'
252    else
253      fileOut = trim(myDir)//trim(listName)//'.results'
254    endif
255    call system('rm -f '//trim(fileOut))  ! clear output file, if it exists
256    do iCase = 1,nCases
257      if (caseType(iCase)=='list') then
258        fileIn = trim(caseDir(iCase))//trim(caseName(iCase))//'.results'
259        call system('cat '//trim(fileIn)//' >> '//trim(fileOut))
260        call system('echo '' '' >> '//trim(fileOut))
261      else
262        call readResult( caseDir(iCase), myRequest, caseName(iCase), &
263                         results=results )
264        open(unitOut,file=trim(fileOut),position='append',status='unknown')
265        write(unitOut,'(10e18.9)') results(1:nResults)
266        close(unitOut)
267      endif
268    end do
269  endif ! (task=='get')
270
271  ! Count jobs and list (but not if it is a list of lists)
272  totJobs = totJobs+nJobs
273  totCores = totCores+nCores
274  if (nJobs>0) totLists = totLists+1
275
276end subroutine scanList
277
278!------------------------------------------------------------------------------
279
280subroutine getCores( queue, ncores )
281
282! Extracts an integer number from string 'queue', assuming that it is the
283! number of cores per job, or returns 1 if there is no such integer
284
285  implicit none
286  character(len=*),intent(in) :: queue    ! queue statement string
287  integer,         intent(out):: ncores   ! number of cores per job
288
289  character(len=1),parameter:: separator(1) = (/' '/)
290  character(len=wl):: words(maxWords)
291  integer:: iWord, n, nWords
292
293  ! Parse queue into blank-separated words
294  call parser(separator,queue,words,nWords)
295
296  ! Find a word that is an integer number. Assume that it is the num. of cores
297  ncores = 1                               ! default value
298  do iWord = 1,nWords
299    n = verify(trim(words(iWord)),'0123456789')  ! verify is an intrinsic funct
300    if (n==0) then                         ! all characters are numbers
301      read(words(iWord),*) ncores          ! read ncores from word
302      return
303    endif
304  end do
305
306end subroutine getCores
307
308!------------------------------------------------------------------------------
309
310subroutine nameJob( jobLine, jobName )
311
312! Generates a job name, concatenating words in a job-specification line
313! .fdf suffixes are removed from words
314
315  implicit none
316  character(len=*),intent(in) :: jobLine   ! string of job specifications
317  character(len=*),intent(out):: jobName   ! job name
318
319  character(len=1),parameter:: separators(2) = (/' ',';'/)
320  integer:: iWord, nc, nWords
321  character(len=wl):: line, word, words(maxWords)
322
323  ! Parse job line
324  line = jobLine
325  call parser(separators,line,words,nWords)
326
327  ! Set job name by concatenating words
328  jobName = ' '
329  do iWord = 1,nWords
330    word = words(iWord)
331    nc = len(trim(word))                         ! number of characters in word
332    if (word(nc-3:nc)=='.fdf') word=word(1:nc-4) ! remove '.fdf' from word
333    jobName = trim(jobName)//'_'//word           ! add word to job name
334  end do ! iWord
335  jobName=jobName(2:)                            ! remove leading '_'
336
337end subroutine nameJob
338
339!------------------------------------------------------------------------------
340
341subroutine runOneJob( dir, queue, files, jobLine )
342
343! Submits (queues) one job, specified by string 'jobLine' of specifications,
344! using the 'queue' statement, within a job-specific subdirectory of 'dir'
345! The subdirectory is named by concatenating the words in 'jobLine' (without
346! '.fdf' suffixes), separated by '_'
347
348  implicit none
349  character(len=*),intent(in) :: dir       ! work directory
350  character(len=*),intent(in) :: queue     ! queuing statement
351  character(len=*),intent(in) :: files     ! data files
352  character(len=*),intent(in) :: jobLine   ! line of job specification
353
354  character(len=1),parameter:: separator(1) = (/';'/)
355  logical:: jobEnded
356  integer:: ic, iLine, iWord, jc, jWord, nc, nLines, nWords
357  character(len=wl):: fileName, jobDir, jobName, line(maxWords), &
358                      myLine, queueJob, word, words(maxWords)
359
360  ! Find job name
361  call nameJob( jobLine, jobName )
362
363  ! Create a new directory for this job and copy files to it
364  jobDir = trim(dir) // trim(jobName) // '/'
365  call system('mkdir -p ' // trim(jobDir))
366  call system('cp -f ' // trim(files) // ' ' // trim(jobDir) // &
367              ' 2> /dev/null')
368
369  ! Parse job-specifications line
370  myLine = jobLine
371  call parser(separator,myLine,words,nWords)
372
373  ! Convert job specifications to fdf format lines
374  nLines = nWords                      ! lines in the output fdf file
375  do iWord = 1,maxWords                ! loop on job specifications
376    word = words(iWord)
377    nc = len(trim(word))
378    if (word(nc-3:nc)=='.fdf') then    ! include new .fdf file
379      line(iWord) = '%include '//trim(word)
380    else                               ! copy statement to new line
381      line(iWord) = word
382    end if
383  end do ! iWord
384
385  ! Write .fdf file for job, but avoid overwriting it (this would occur only
386  ! if job contains a single .fdf file, with no other specifications)
387  fileName = trim(jobName)//'.fdf'
388  if (trim(fileName)/=trim(words(1))) then
389    open(unitOut,file=trim(jobDir)//trim(fileName))
390    write(unitOut,'(a)') 'SystemLabel '//trim(jobName) ! this is siesta-specific
391    do iLine = nLines,1,-1               ! last lines (words) have priority
392      write(unitOut,'(a)') trim(line(iLine))
393    end do
394    close(unitOut)
395  endif
396
397  ! Set job
398  queueJob = queue
399  nc = len(trim(queueJob))
400  jc = len('$jobName')
401  do ic = nc-jc+1,1,-1
402    if (queueJob(ic:ic+jc-1)=='$jobName') &
403      queueJob(ic:) = trim(jobName)//queueJob(ic+jc:)
404  end do
405!  print*, 'runOneJob: queueJob = ',trim(queueJob)
406
407  ! Submit job only if directory does not contain terminated results already
408  ! This is intended to re-run a whole list for failed jobs
409
410  call chdir(trim(jobDir))
411
412  ! Check for existence of a special file which signals the end of execution
413  inquire(file=trim(file0Exit),exist=jobEnded)
414
415  if (.not.jobEnded) then
416     ! Fall back to old test with EIG file, in case of an older
417     ! version of Siesta without the new flag file.
418     inquire(file=trim(jobName)//endSuffix,exist=jobEnded)
419  endif
420
421  if (.not.jobEnded) then
422     !print "(2a)", 'queueJob = ',trim(queueJob)
423     call system(trim(queueJob))
424  endif
425  call chdir(trim(dir))
426
427end subroutine runOneJob
428
429!------------------------------------------------------------------------------
430
431subroutine readLine( unit, line, iostat )
432
433! Reads one line of a job list file (specified by a I/O unit), removing
434! comments (marked by a leading '#'), and including continuation lines
435! (marked by a trailing '\')
436
437  implicit none
438  integer,         intent(in) :: unit    ! input file unit
439  character(len=*),intent(out):: line    ! line read
440  integer,optional,intent(out):: iostat  ! I/O status
441
442  logical keepReading
443  integer:: lc, nc, newc, status
444  character:: lastChar
445  character(len=1024):: newLine    ! a sufficiently long array
446
447  line = ' '                       ! set a blank line to start
448  nc = 0                           ! number of characters in line
449  keepReading = .true.
450  do while (keepReading)           ! loop on continuation lines
451
452    ! Read new line from input file
453    read(unit,'(a)',iostat=status) newLine
454
455    ! End-of-file trap
456    if (present(iostat)) iostat=status
457    if (status<0) then
458      line = ' '
459      return
460    endif
461
462    ! Remove comments and leading blanks
463    lc = scan(newLine,'#')-1             ! last character before comments
464    if (lc>=0) newLine = newLine(1:lc)   ! select characters before comments
465    newLine = adjustl(newLine)           ! remove leading blanks
466
467    ! Find last character and remove continuation mark
468    newc = len(trim(newLine))            ! nonblank characters in new line
469    lastChar = newLine(newc:newc)        ! last nonblank character
470    if (lastChar==char(92)) then         ! line will continue  (backslash) '\'
471      newLine(newc:newc) = ' '           ! remove '\' mark
472      newc = len(trim(newLine))          ! remaining nonblank characters
473      keepReading = .true.
474    else
475      keepReading = .false.
476    endif
477
478    ! Concatenate continued line
479    if (nc==0) then                      ! this is first line
480      line = newLine(1:newc)
481      nc = newc
482    else                                 ! this is a continuation line
483      line = line(1:nc)//' '//newLine(1:newc)
484      nc = nc+1+newc
485    endif
486
487    ! Check array length
488    if (nc>len(line)) stop 'readLine ERROR: len(line) too small'
489
490  end do ! while (keepReading)
491
492end subroutine readLine
493
494!------------------------------------------------------------------------------
495
496subroutine parser( separators, line, words, nWords )
497
498! Parses one line into words, using any of a number of separators.
499! If ' ' is not one of the separators, the 'words' may contain blanks,
500! but leading and trailing blanks will be removed from them in any case
501
502  implicit none
503  character(len=*),intent(in)   :: separators(:) ! character(s) between words
504  character(len=*),intent(in)   :: line      ! line to be parsed
505  character(len=*),intent(out)  :: words(:)  ! words in line, without comments
506  integer,         intent(out)  :: nWords    ! number of words
507
508  character(len=ll):: myLine
509  integer:: is, iWord, nc, ncs, ns
510
511  ! Loop on words
512  myLine = adjustl(line)                  ! copy line, removing leading blanks
513  ns = size(separators)                   ! number of alternative separators
514  do iWord = 1,size(words)                ! avoid overflooding array size
515    nc = len(trim(myLine))                ! number of nonblank characters
516    do is = 1,ns                          ! loop on separators
517      ncs = scan(myLine,separators(is))-1 ! last character before this separator
518      if (ncs>0) nc = min(nc,ncs)         ! character before first separator
519    end do
520    if (nc<=0) then                       ! no more words
521      nWords = iWord-1
522      return                              ! normal return point
523    else
524      if (nc>len(words(iWord))) stop 'parser ERROR: len(words) too small'
525      words(iWord) = adjustl(myLine(1:nc))
526      myLine = adjustl(myLine(nc+2:))     ! discard previous word from myLine
527    endif
528  end do ! iWord
529  stop 'parser ERROR: size(words) too small'
530
531end subroutine parser
532
533!------------------------------------------------------------------------------
534
535subroutine readResult( dir, request, name, result, results )
536
537! Reads magnitudes specified by 'request' string(s), calculated by a SIESTA
538! job named 'name', within directory 'dir'. The magnitudes are output in one
539! or both of: derived-type structure 'result'; and real array 'results'
540! Presently allowed values (case sensitive) of 'request' are
541!   (volume|energy|pressure|virial|maxForce|avgForce)
542! Additionally, derived-type 'result' returns the system geometry, forces,
543! and stress tensor
544
545  implicit none
546  character(len=*),          intent(in) :: dir         ! job directory
547  character(len=*),          intent(in) :: request(:)  ! requested resuts
548  character(len=*),          intent(in) :: name        ! job name
549  type(resultsType),optional,intent(out):: result      ! job-results structure
550  real(dp),         optional,intent(out):: results(:)  ! requested job results
551
552  integer :: ia, ic, is, iostat, iResult, nAtoms, nResults, za
553  real(dp):: c(3,3)
554  character(len=wl):: fileName
555  type(resultsType) :: res
556
557  ! Read final structure
558  fileName = trim(dir)//trim(name)//'.XV'
559  open(unitIn,file=trim(fileName),status='old',iostat=iostat)
560!  print'(a,i3,2x,a)','readResult: iostat,fileName=', iostat, trim(fileName)
561  if (iostat==0) then
562    do ic = 1,3
563      read(unitIn,*) res%cell(:,ic)
564    end do
565    read(unitIn,*) nAtoms
566    allocate( res%za(nAtoms), res%coord(3,nAtoms) )
567    do ia = 1,nAtoms
568      read(unitIn,*) is, res%za(ia), res%coord(:,ia)
569    end do
570  else
571    nAtoms = 0
572    allocate( res%za(nAtoms), res%coord(3,nAtoms) )
573    res%cell = 0
574    res%coord = 0
575    res%za = 0
576  end if
577  close(unitIn)
578
579  ! Read energy, force, and stress
580  fileName = trim(dir)//'FORCE_STRESS'
581  open(unitIn,file=trim(fileName),status='old',iostat=iostat)
582  if (iostat==0) then
583    read(unitIn,*) res%energy
584    read(unitIn,*) res%stress
585    read(unitIn,*) nAtoms
586    deallocate(res%za)
587    allocate(res%force(3,nAtoms),res%za(nAtoms))
588    do ia = 1,nAtoms
589      read(unitIn,*) is, res%za(ia), res%force(:,ia)
590    end do
591  else
592    allocate( res%force(3,0) )
593    res%energy = 0
594    res%stress = 0
595    res%force = 0
596  end if
597  close(unitIn)
598
599  ! Find cell volume, pressure, and virial
600  c = res%cell
601  res%volume = abs( c(1,1)*c(2,2)*c(3,3) - c(1,1)*c(2,3)*c(3,2) + &
602                    c(1,2)*c(2,3)*c(3,1) - c(1,2)*c(2,1)*c(3,3) + &
603                    c(1,3)*c(2,1)*c(3,2) - c(1,3)*c(2,2)*c(3,1) )
604  res%pressure = -res%stress(1,1)-res%stress(2,2)-res%stress(3,3)
605  res%virial = sum(res%coord*res%force)
606
607  ! Copy results to output structure
608  if (present(result)) then
609    result = res
610  end if
611
612  ! Select requested (scalar) results
613  if (present(results)) then
614    nResults = count(request/=' ')
615    do iResult = 1,nResults
616      select case (trim(request(iResult)))
617        case ('energy')
618          results(iResult) = res%energy
619        case ('volume')
620          results(iResult) = res%volume
621        case ('pressure')
622          results(iResult) = res%pressure
623        case ('virial')
624          results(iResult) = res%virial
625        case ('maxForce')
626          results(iResult) = sqrt(maxval(sum(res%force**2,1)))
627        case ('avgForce')
628          results(iResult) = sqrt(sum(res%force**2)/max(1,nAtoms))
629      end select
630    end do
631  end if
632
633  ! Deallocate internal pointers
634  if (associated(res%za)) deallocate(res%za)
635  if (associated(res%coord)) deallocate(res%coord)
636  if (associated(res%force)) deallocate(res%force)
637
638end subroutine readResult
639
640END MODULE jobList
641