text:
      SUBROUTINE STRM61
C
C
C   PHASE I OF  STRESS DATA RECOVERY FOR TRIANGULAR MEMBRANE ELEMENT TRI
C
C   OUTPUTS FROM THIS PHASE FOR USE IN PHASE II ARE THE FOLLOWING
C
C   1) ELEMENT ID             WORDS    1    STORAGE IN PH1OUT   1
C   2) SIX SILS               WORDS    6                        2-7
C   3) THICKNESS T1           WORDS    1                        8
C   4) THICKNESS T2           WORDS    1                        9
C   5) THICKNESS T3           WORDS    1                       10
C   6) REFERENCE TEMP T0      WORDS    1                       11
C   7) S SUB I MATRICES       WORDS    216                     12-227
C   8) THERMAL VECTOR G ALF   WORDS    3                      228-230
C
C    EST ENTRIES SAME AS IN SUBROUTINE KTRM6S
C
C
      REAL NSM,IVECT,JVECT,KVECT
C
      DIMENSION IEST(45),IND(6,3),EE1(6),NPH1OU(990),XC(6),YC(6),ZC(6)
     1,   Q(6,6),QQ(36),IVECT(3),JVECT(3),KVECT(3),E( 6 ),EPH1(6)
     2,   NAME(2),ICS(6),NL(6),TRANS(9),BALOTR(9),EMOD(9),TM(3,12)
     3,   TMM(36)
C
      COMMON /SDR2X5/ EST(100),PH1OUT(250)
      COMMON /MATIN / MATID,MATFLG,ELTEMP,PLA34,SINTH,COSTH
      COMMON /MATOUT/ EM(6),RHOY,ALF(3),TREF,GSUBE,SIGTY,SIGCY,SIGSY
C
      EQUIVALENCE (NPH1OU(1),PH1OUT(1)),(IEST(1),EST(1))
      EQUIVALENCE (TM(1,1),TMM(1))
C
      DATA NAME /4HSTRM,4H61  / , BLANK  /4H    /
      DATA DEGRA /0.0174532925/
C
      IDELE=IEST(1)
      DO 109 I=1,6
      NL(I)=IEST(I+1)
  109 CONTINUE
      THETAM=EST(8)
      MATID1=IEST(9)
      TMEM1 =EST(10)
      TMEM3 =EST(11)
      TMEM5 =EST(12)
C
C   IF TMEM3 OR TMEM5 IS 0.0 OR BLANK , IT WILL BE SET EQUAL TO TMEM1
C
      IF (TMEM3.EQ.0.0. OR .TMEM3.EQ.BLANK)  TMEM3 = TMEM1
C
      IF (TMEM5.EQ.0.0. OR .TMEM5.EQ.BLANK)  TMEM5 = TMEM1
C
      NSM = EST(13)
C
      J=0
      DO 120 I=14,34,4
      J=J+1
      ICS(J)=IEST(I)
      XC(J) = EST(I+1)
      YC(J) = EST(I+2)
      ZC(J)=EST(I+3)
  120 CONTINUE
      ELTEMP=(EST(38)+EST(39)+EST(40)+EST(41)+EST(42)+EST(43))/6.0
      THETA1=THETAM*DEGRA
      SINTH=SIN(THETA1)
      COSTH=COS(THETA1)
      IF (ABS(SINTH).LE.1.0E-06) SINTH=0.0
C
C
C   EVALUATE MATERIAL PROPERTIES
C
      MATFLG = 2
      MATID  = MATID1
      CALL MAT (IDELE)
      TO=TREF
C
C   CALCULATIONS FOR THE TRIANGLE
C
      CALL TRIF (XC,YC,ZC,IVECT,JVECT,KVECT,DISTA,DISTB,DISTC,IEST(1),
     1           NAME)
C
C   TRANSFORMATION MATRIX BETWEEN ELEMENT AND BASIC CO-ORDINATES
C
      E(1)=IVECT(1)
      E(2)=JVECT(1)
      E(3)=IVECT(2)
      E(4)=JVECT(2)
      E(5)=IVECT(3)
      E(6)=JVECT(3)
C
C   CALCULATIONS FOR  Q MATRIX AND ITS INVERSE
C
      DO 110 I=1,6
      DO 110 J=1,6
      Q(I,J)=0.0
  110 CONTINUE
      DO 115 I=1,6
      Q(I,1)=1.0
      Q(I,2)=XC(I)
      Q(I,3)=YC(I)
      Q(I,4)=XC(I)*XC(I)
      Q(I,5)=XC(I)*YC(I)
      Q(I,6)=YC(I)*YC(I)
  115 CONTINUE
C
C     FIND INVERSE OF Q  MATRIX
C
C     NO NEED TO COMPUTE DETERMINANT SINCE IT IS NOT USED SUBSEQUENTLY.
      ISING = -1
      CALL INVERS (6,Q,6,QQ(1),0,DETERM,ISING,IND)
C
C   ISING EQUAL TO 2 IMPLIES THAT Q MATRIX IS SINGULAR
C
      DO 152 I=1,6
      DO 152 J=1,6
      IJ=(I-1)*6+J
      QQ(IJ)=Q(I,J)
  152 CONTINUE
      DO 154 I=1,9
      BALOTR(I)=0.0
  154 CONTINUE
C
      DO 102 I=1,7
      PH1OUT(I)=EST(I)
  102 CONTINUE
      PH1OUT(8)=EST(10)
      PH1OUT(9)=EST(11)
      PH1OUT(10)=EST(12)
      PH1OUT(11)=TO
      EMOD(1)=EM(1)
      EMOD(2)=EM(2)
      EMOD(3)=EM(3)
      EMOD(4)=EM(2)
      EMOD(5)=EM(4)
      EMOD(6)=EM(5)
      EMOD(7)=EM(3)
      EMOD(8)=EM(5)
      EMOD(9)=EM(6)
C
C   STRESSES AND STRAINS ARE EVALUATED AT FOUR POINTS ,VIZ., THE THREE
C   CORNER GRID POINTS AND THE CENTROID
C
      DO 700 JJ=1,4
      J=2*(JJ-1)+1
      IF (J.EQ.7) GO TO 103
      X=XC(J)
      Y=YC(J)
      GO TO 104
  103 X=(XC(1)+XC(3)+XC(5))/3.0
      Y=(YC(1)+YC(3)+YC(5))/3.0
  104 CONTINUE
      DO 105 I=1,36
      TMM(I)=0.0E0
  105 CONTINUE
C
C   TM MATRIX IS THE PRODUCT OF B AND QINVERSE MATRICES
C
      DO 258 J=1,6
      J1=(J-1)*2+1
      J2=J1+1
      TM(1,J1)=Q(2,J)+2.0*X*Q(4,J)+Y*Q(5,J)
      TM(2,J2)=Q(3,J)+X*Q(5,J)+2.0*Y*Q(6,J)
      TM(3,J1)=TM(2,J2)
      TM(3,J2)=TM(1,J1)
  258 CONTINUE
      DO 600 II=1,6
      IF (ICS(II).EQ.0) GO TO 130
      CALL TRANSS (IEST(4*II+10),TRANS)
      CALL GMMATS (E,3,2,+1,TRANS,3,3,0,EE1)
      GO TO 133
  130 CONTINUE
      DO 132 I=1,3
      DO 132 J=1,2
      I1=(I-1)*2+J
      J1=(J-1)*3+I
      EE1(J1)=E(I1)
  132 CONTINUE
  133 CONTINUE
      IJ1=(JJ-1)*54+(II-1)*9+12
      MZ=(II-1)*6+1
      CALL GMMATS (EMOD,3,3,0,TMM(MZ),2,3,+1,EPH1)
      CALL GMMATS (EPH1,3,2,0,EE1,2,3,0,PH1OUT(IJ1))
  600 CONTINUE
  700 CONTINUE
      CALL GMMATS (EMOD,3,3,0,ALF,3,1,0,PH1OUT(228))
      RETURN
      END


text:
c
c --------------------------------------------------------------
c
      subroutine preintcopy(val,mitot,nvar,ilo,ihi,
     1                      level,fliparray)
c
      use amr_module
      implicit double precision (a-h,o-z)

      dimension fliparray((mitot)*nghost*nvar)
      dimension val(nvar,mitot)
      dimension ist(3), iend(3), ishift(3)

c  NEW INDEXING ORDER SWITCHED
      iadd(ivar,i)  = locflip + ivar-1 + nvar*(i-1)
c
c  :::::::::::::: PREINTCOPY :::::::::::::::::::::::::::::::::::::::::::
c     For periodic boundary conditions more work needed to initialize a
c     new grid that sticks out. This routine was
c     called because the patch sticks out of the domain,
c     and has periodic bc.s preprocess the patch before calling
c     intcopy to shift the patch periodically back into the domain.
c
c     Inputs to this routine:
c     ilo,ihi = the location in index space of this patch.
c
c     Outputs from this routine:
c     The values of the grid are inserted
c     directly into the enlarged val array for this piece.
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
        locflip = 1

c
c     # will divide patch into 3 possibilities (some empty):
c       x sticks out left, x interior, x sticks out right
c       for example, the max. would be
c       i from (ilo,-1), (0,iregsz(level)-1), (iregsz(level),ihi)

        if (xperdom) then
           ist(1) = ilo
           ist(2) = 0
           ist(3) = iregsz(level)
           iend(1) = -1
           iend(2) = iregsz(level)-1
           iend(3) = ihi
           ishift(1) = iregsz(level)
           ishift(2) = 0
           ishift(3) = -iregsz(level)
        else    ! if not periodic, set vals to only have nonnull intersection for interior region
           ist(1)    = iregsz(level)
           ist(2)    = ilo
           ist(3)    = iregsz(level)
           iend(1)   = -1
           iend(2)   = ihi
           iend(3)   = -1
           ishift(1) = 0
           ishift(2) = 0
           ishift(3) = 0
        endif

       do 20 i = 1, 3
          i1 = max(ilo,  ist(i))
          i2 = min(ihi, iend(i))
          if (i1 .le. i2) then ! part of patch in this region
c
            iputst = (i1 - ilo) + 1
            call intcopy(val,mitot,nvar,
     2                    i1+ishift(i),i2+ishift(i),
     3                    level,iputst)

          endif

 20    continue

      return
      end

ext: f

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


text:
module ioModule
 use varModule

 private

 interface matprint
    module procedure print_eigenproblem
    module procedure print_vector
    module procedure print_vector_with_header
    module procedure print_matrix
    module procedure print_matrix_with_header
 end interface

 interface print_gamess_vecs
    module procedure print_gamess_vec
 end interface

 interface printHeader
     module procedure printHeader
 end interface

 public :: matprint, print_gamess_vecs, printHeader

 contains

 subroutine printHeader(header)
  character(*), intent(in) :: header
  character(len=30)        :: headerfmt

  lspace = 7 + int((50-len(header))/2)
  write(headerfmt,'("(",i2,"x,a)")') lspace
  write(*,*)
  write(*,headerfmt) repeat('=',len(trim(header)))
  write(*,headerfmt) trim(header)
  write(*,headerfmt) repeat('=',len(trim(header)))
 end subroutine printHeader

 subroutine print_eigenproblem(M, V)
  integer(ik), parameter :: vec_per_page = 5
  integer(ik)          :: i, j, iel, k, l
  real(dp), intent(in) :: M(:,:)
  real(dp), intent(in) :: V(:)

 if (mod(size(M,2),vec_per_page) == 0 ) then
   iel = size(M,2)/vec_per_page
 else
   iel = size(M,2)/vec_per_page+1
 endif

 do i = 1, iel
   if(i < iel) then
     write(*,'(/10x,5(i5,5x))') (j+(i-1)*5,j=1,5)
     write(*,'(/7x,5(f10.5)/)') (V(j+(i-1)*10),j=1,5)
       do k = 1,size(M,1)
         write(*,'(1x,i3,3x,5f10.5)') k, (M(k,l+(i-1)*5),l=1,5)
       enddo
       write(*,*)
   else
     write(*,'(/10x,5(i5,5x))') (j+(i-1)*5,j=1,5-(iel*5-size(M,2)))
     write(*,'(/7x,5(f10.5)/)') (V(j+(i-1)*5),j=1,5-(iel*5-size(M,2)))
       do k = 1,size(M,1)
         write(*,'(1x,i3,3x,5f10.5)') k, (M(k,l+(i-1)*5),l=1,5-(iel*5-size(M,2)))
       enddo
   endif
 enddo
 write(*,*)
 end subroutine print_eigenproblem

 subroutine print_vector(V)
  real(dp), intent(in) :: V(:)
  integer(ik)          :: i

  do i = 1, size(V)
    write(*,'(1x,i3,f10.5)') i, V(i)
  enddo
 end subroutine print_vector

 subroutine print_vector_with_header(V, header)
  real(dp),     intent(in) :: V(:)
  character(*), intent(in) :: header
  integer(ik)              :: i

  write(*,'(/4x,a)') repeat('=',len(trim(header)))
  write(*,'(4x,a)') trim(header)
  write(*,'(4x,a/)') repeat('=',len(trim(header)))

  do i = 1, size(V)
    write(*,'(1x,i3,f10.5)') i, V(i)
  enddo
 end subroutine print_vector_with_header

 subroutine print_matrix(M)
  real(dp), intent(in) :: M(:,:)
  integer(ik)          :: i,j,iel,k,l

 if (mod(size(M,2),5) == 0 ) then
   iel = size(M,2)/5
 else
   iel = size(M,2)/5+1
 endif

 do i = 1, iel
   if(i < iel) then
     write(*,'(/10x,5(i5,5x))') (j+(i-1)*5,j=1,5)
       do k = 1,size(M,1)
         write(*,'(1x,i3,3x,5f10.5)') k, (M(k,l+(i-1)*5),l=1,5)
       enddo
       write(*,*)
   else
     write(*,'(/10x,5(i5,5x))') (j+(i-1)*5,j=1,5-(iel*5-size(M,2)))
       do k = 1,size(M,1)
         write(*,'(1x,i3,3x,5f10.5)') k, (M(k,l+(i-1)*5),l=1,5-(iel*5-size(M,2)))
       enddo
   endif
 enddo
 write(*,*)
 end subroutine print_matrix

 subroutine print_matrix_with_header(M, header)
  real(dp),     intent(in) :: M(:,:)
  character(*), intent(in) :: header
  integer(ik)              :: i,j,iel,k,l
  character(len=30)        :: headerfmt

  lspace = 7 + int((50-len(header))/2)
  write(headerfmt,'("(",i2,"x,a)")') lspace
  write(*,headerfmt) repeat('=',len(trim(header)))
  write(*,headerfmt) trim(header)
  write(*,headerfmt) repeat('=',len(trim(header)))

 if (mod(size(M,2),5) == 0 ) then
   iel = size(M,2)/5
 else
   iel = size(M,2)/5+1
 endif

 do i = 1, iel
   if(i < iel) then
     write(*,'(/10x,5(i5,5x))') (j+(i-1)*5,j=1,5)
       do k = 1,size(M,1)
         write(*,'(1x,i3,3x,5f10.5)') k, (M(k,l+(i-1)*5),l=1,5)
       enddo
       write(*,*)
   else
     write(*,'(/10x,5(i5,5x))') (j+(i-1)*5,j=1,5-(iel*5-size(M,2)))
       do k = 1,size(M,1)
         write(*,'(1x,i3,3x,5f10.5)') k, (M(k,l+(i-1)*5),l=1,5-(iel*5-size(M,2)))
       enddo
   endif
 enddo
 write(*,*)
 end subroutine print_matrix_with_header

 subroutine print_gamess_vec(coeffs)
  real(dp),           intent(in) :: coeffs(:,:)
  !character(len=100), intent(in) :: vec_file
  integer(ik) :: i, j, l, ilab, llab, nlines
  integer(ik), parameter :: gunit=16

  if (mod(size(coeffs,1),5) == 0) then
    nlines = int(size(coeffs,1)/5)
  else
    nlines = int(size(coeffs,1)/5,ik)+1
  endif

  print *, 'nlines = ', nlines

  open(gunit, file='vectors.out', status="replace")
 ! write header and then orbitals
  write(gunit,'(/a,i3,a/)') ' $guess guess=moread norb=',size(coeffs,2),' punmo=.true. prtmo=.true. $end'
  write(gunit,'(a)') ' $vec'
  do i = 1, size(coeffs,2)
    if (i >= 100) then
      ilab = mod(i, 100)
    else
      ilab = i
    endif
    do l = 1, nlines
      if (l > 1000) then
        llab = mod(l, 1000)
      else
        llab = l
      endif
      if (l < nlines) then
        write(gunit,gmsvfmt) ilab,llab,(coeffs(j+(l-1)*5,i),j=1,5)
      else
        write(gunit,gmsvfmt) ilab,llab,(coeffs(j+(l-1)*5,i),j=1,5-(l*5-size(coeffs,1)))
      endif
    enddo
  enddo
  write(gunit,'(a)') ' $end'
  close(gunit)

 end subroutine print_gamess_vec

end module ioModule

max_stars_repo_name: lmmentel/gamess-us-fortran-interface

max_issues_repo_name: lmmentel/gamess-us-fortran-interface

max_forks_repo_name: lmmentel/gamess-us-fortran-interface


text:
#include "error.inc"
program vasp_driver
use libatoms_module
implicit none

   type(Atoms) :: at

   character(len=STRING_LENGTH) :: infile, outfile
   character(len=STRING_LENGTH) :: args_str
   character(len=STRING_LENGTH)  :: arg
   integer :: i, index_insert
   type(CInOutput) :: outfile_io

   integer :: error, stat_vals(13)

   call system_initialise(verbosity=PRINT_SILENT, enable_timing=.true.)
   call verbosity_push(PRINT_NORMAL)
   call system_timer('vasp_filepot')

   if (cmd_arg_count() < 2) then
      call system_abort("Usage: vasp_filepot infile outfile [ other_cli_arg1=v1 other_cli_arg2=v2 ...]")
   endif

   call get_cmd_arg(1, infile)
   call get_cmd_arg(2, outfile)

    args_str = ""
    if (cmd_arg_count() > 2) then
      do i=3, cmd_arg_count()
        call get_cmd_arg(i, arg)
        !add {} if there is space in the arg
        if (index(trim(arg)," ").ne.0) then
            index_insert = index(trim(arg),"=")
            arg(index_insert+1:len_trim(arg)+2) = "{"//arg(index_insert+1:len_trim(arg))//"}"
            call print('arg: '//trim(arg),PRINT_SILENT)
        endif
        args_str = trim(args_str) // " " // trim(arg)
      end do
    endif

    call read(at,infile)

    call print('vasp_filepot args_str '//trim(args_str), PRINT_ALWAYS)

    call do_vasp_calc(at, trim(args_str), error)
    HANDLE_ERROR(error)

    call initialise(outfile_io,outfile,action=OUTPUT)
    call write(at, outfile_io, real_format='%20.13f')
    call finalise(outfile_io)

    call finalise(at)

    call system_timer('vasp_filepot')
    mainlog%prefix=""
    call verbosity_push(PRINT_SILENT)
    call system_finalise()

contains

subroutine do_vasp_calc(at, args_str, error)
   type(Atoms), intent(inout) :: at
   character(len=*), intent(in) :: args_str
   integer, intent(out), optional :: error

   character(len=STRING_LENGTH) :: incar_template_file, kpoints_file, potcar_files, vasp_path, run_suffix, verbosity_str
   character(len=STRING_LENGTH) :: calc_energy, calc_force, calc_virial, calc_local_energy
   logical :: do_calc_energy, do_calc_force, do_calc_virial, do_calc_local_energy, persistent, persistent_n_atoms_exist
   logical :: clean_up_files, ignore_convergence, no_use_WAVECAR, force_constant_basis
   integer :: clean_up_keep_n
   integer :: persistent_restart_interval, persistent_n_atoms
   real(dp) :: persistent_max_wait_time
   type(Dictionary) :: incar_dict, cli

   integer :: persistent_run_i
   logical :: persistent_already_started = .false.
   type(inoutput) :: persistent_io

   integer :: run_dir_i, force_run_dir_i

   integer :: i_charge
   real(dp) :: nelect, nelect_of_elem(size(ElementName)), r_charge
   real(dp), pointer :: magmom_p(:)

   character(len=STRING_LENGTH) :: run_dir
   type(Inoutput) :: io

   integer :: at_i, stat
   integer, pointer :: sort_index_p(:)
   logical :: converged, have_wavecar

   INIT_ERROR(error)

   call initialise(cli)
   call param_register(cli, "verbosity", "NORMAL", verbosity_str, help_string="verbosity level")
   call param_register(cli, "INCAR_template", "INCAR.template", incar_template_file, help_string="name of file with template for main input (INCAR) file")
   call param_register(cli, "kpoints_file", "KPOINTS", kpoints_file, help_string="name of k-points file (KPOINTS)")
   call param_register(cli, "potcar_files", PARAM_MANDATORY, potcar_files, help_string="mapping from atomic numbers to POTCAR files, Z_1 POTCAR_FILE_1 [ Z_2 POTCAR_FILE_2 ... ]")
   call param_register(cli, "vasp", PARAM_MANDATORY, vasp_path, help_string="vasp executable")
   call param_register(cli, "energy", "", calc_energy, help_string="if set, calculate energy and put it in parameter named this")
   call param_register(cli, "force", "", calc_force, help_string="if set, calculate forces and put them in property named this")
   call param_register(cli, "virial", "", calc_virial, help_string="if set, calculate virial and put it in parameter named this")
   call param_register(cli, "local_energy", "", calc_local_energy, help_string="not supported")
   call param_register(cli, "clean_up_files", "T", clean_up_files, help_string="if true, delete run directory when done")
   call param_register(cli, "clean_up_keep_n", "1", clean_up_keep_n, help_string="if cleaning up, keep this many directories around")
   call param_register(cli, "ignore_convergence", "F", ignore_convergence, help_string="If true, accept result even if failed to converge SCF")
   call param_register(cli, "no_use_WAVECAR", "F", no_use_WAVECAR, help_string="If true, don't reuse previous wavefunction file WAVECAR")
   call param_register(cli, "force_constant_basis", "F", force_constant_basis, help_string="If true, do a restart with the same basis")
   call param_register(cli, "run_suffix", "run", run_suffix, help_string="suffix to label persistent state from run to run")
   call param_register(cli, "persistent", "F", persistent, help_string="If true, use persistent connection to a long-running vasp process")
   call param_register(cli, "persistent_restart_interval", "100", persistent_restart_interval, help_string="If persistent, restart vasp this often")
   call param_register(cli, "persistent_max_wait_time", "600", persistent_max_wait_time, help_string="If persistent, never wait more than this long (s) for calculation to finish")
   if (.not. param_read_line(cli, args_str, ignore_unknown=.true.,task='do_vasp_calc args_str')) then
      call print("Args:  verbosity INCAR_template kpoints_file potcar_files vasp energy force", PRINT_ALWAYS)
      call print("       virial local_energy clean_up_files ignore_convergence no_use_WAVECAR", PRINT_ALWAYS)
      call print("       force_constant_basis run_suffix", PRINT_ALWAYS)
      RAISE_ERROR('do_vasp_calc failed to parse args_str="'//trim(args_str)//'"', error)
   endif
   call finalise(cli)

   call verbosity_push(verbosity_of_str(trim(verbosity_str)))

   do_calc_energy = len_trim(calc_energy) > 0
   do_calc_force = len_trim(calc_force) > 0
   do_calc_virial = len_trim(calc_virial) > 0
   do_calc_local_energy = len_trim(calc_local_energy) > 0

   call print("do_vasp_calc using args:", PRINT_VERBOSE)
   call print("  INCAR_template " // trim(INCAR_template_file), PRINT_VERBOSE)
   call print("  kpoints_file " // trim(kpoints_file), PRINT_VERBOSE)
   call print("  potcar_files " // trim(potcar_files), PRINT_VERBOSE)
   call print("  vasp " // trim(vasp_path), PRINT_VERBOSE)
   call print("  energy " // do_calc_energy, PRINT_VERBOSE)
   call print("  force " // do_calc_force, PRINT_VERBOSE)
   call print("  virial " // do_calc_virial, PRINT_VERBOSE)
   call print("  local_energy " // do_calc_local_energy, PRINT_VERBOSE)
   call print("  clean_up_files " // clean_up_files, PRINT_VERBOSE)
   call print("  clean_up_keep_n " // clean_up_keep_n, PRINT_VERBOSE)
   call print("  ignore_convergence " // ignore_convergence, PRINT_VERBOSE)
   call print("  no_use_WAVECAR " // no_use_WAVECAR, PRINT_VERBOSE)
   call print("  force_constant_basis " // force_constant_basis, PRINT_VERBOSE)
   call print("  run_suffix " // run_suffix, PRINT_VERBOSE)
   call print("  persistent " // persistent, PRINT_VERBOSE)
   if (persistent) then
      call print("  persistent_restart_interval " // persistent_restart_interval, PRINT_VERBOSE)
      call print("  persistent_max_wait_time " // persistent_max_wait_time, PRINT_VERBOSE)
   endif

   call system_timer('vasp_filepot/prep')
   ! check for local energy
   if (do_calc_local_energy) then
      RAISE_ERROR("do_vasp_calc can't do local_energy", error)
   endif

   call print("reading template", PRINT_VERBOSE)
   ! read the template incar
   call read_vasp_incar_dict(incar_dict, trim(incar_template_file), error)
   PASS_ERROR(error)

   call print("setting INCAR dictionary parameters", PRINT_VERBOSE)
   ! check for stress, set ISIF
   if (.not. do_calc_energy .and. .not. do_calc_force .and. .not. do_calc_virial) then
      return
   else if (.not. do_calc_virial) then
      call set_value(incar_dict, "ISIF", 0)
   else
      if (persistent) then
         call set_value(incar_dict, "ISIF", 3)
      else
         call set_value(incar_dict, "ISIF", 2)
      endif
   endif

   ! remove EDIFFG to avoid VASP quitting when it thinks convergence has been reached
   if (has_key(incar_dict,"EDIFFG")) call remove_value(incar_dict,"EDIFFG")

   force_run_dir_i = -1
   if (persistent) then ! make run_dir manually
      inquire(file="persistent_run_i", exist=persistent_already_started)
      run_dir = "vasp_run_0"
      force_run_dir_i = 0
      if (.not. persistent_already_started) then
         call system_command("mkdir "//trim(run_dir))
         call initialise(persistent_io, "persistent_run_i", action=OUTPUT)
         call print("-1", file=persistent_io)
         call finalise(persistent_io)
      endif

      inquire(file="persistent_n_atoms", exist=persistent_n_atoms_exist)
      if (persistent_n_atoms_exist) then
         call initialise(persistent_io, "persistent_n_atoms", action=INPUT)
         line = read_line(persistent_io)
         read (unit=line, fmt=*, iostat=stat) persistent_n_atoms
         call finalise(persistent_io)
      else
         persistent_n_atoms = 0
      end if
      call print('got old persistent_n_atoms = '//persistent_n_atoms, PRINT_VERBOSE)
      call print('writing new persistent_n_atoms = '//at%n, PRINT_VERBOSE)
      call initialise(persistent_io, "persistent_n_atoms", action=OUTPUT)
      call print(at%n, file=persistent_io)
      call finalise(persistent_io)
   endif
   run_dir = make_run_directory("vasp_run", force_run_dir_i, run_dir_i)

   if (persistent) then
   ! set dynamics to REFTRAJ, number of ionic steps to 1000000
      call set_value(incar_dict, "NSW", 1000000)
      call set_value(incar_dict, "IBRION", 12)
   else
   ! set dynamics to MD, number of ionic steps to 0
      call set_value(incar_dict, "NSW", 0)
      call set_value(incar_dict, "IBRION", 0)
   endif
   if (no_use_WAVECAR) then
      call set_value(incar_dict, "ISTART", 0)
   else
      inquire(file="WAVECAR."//trim(run_suffix), exist=have_wavecar)
      if (have_wavecar) then
         call print("copying old WAVECAR.run_suffix to run_dir/WAVECAR", PRINT_VERBOSE)
         call system_command("cp WAVECAR."//trim(run_suffix)//" "//trim(run_dir)//"/WAVECAR", status=stat)
         if (stat /= 0) then
            RAISE_ERROR("do_vasp_calc copying file WAVECAR."//trim(run_suffix)//" into "//trim(run_dir), error)
         endif
      endif
      if (force_constant_basis) then
         ! 3 is faster, but not tested
         call set_value(incar_dict, "ISTART", 2)
      else
         call set_value(incar_dict, "ISTART", 1)
      endif
   endif

   call print("doing sorting", PRINT_VERBOSE)
   call add_property(at, "sort_index", 0, n_cols=1, ptr=sort_index_p, error=error)
   PASS_ERROR(error)
   do at_i=1, at%N
      sort_index_p(at_i) = at_i
   end do
   call atoms_sort(at, 'Z', error=error)
   PASS_ERROR(error)

   if (persistent) then
      call print("reading persistent_run_i", PRINT_VERBOSE)
      call initialise(persistent_io, "persistent_run_i", action=INPUT)
      line = read_line(persistent_io)
      read (unit=line, fmt=*, iostat=stat) persistent_run_i
      call finalise(persistent_io)
      call print("got persistent_run_i="//persistent_run_i, PRINT_VERBOSE)
      persistent_run_i = persistent_run_i + 1
      call initialise(persistent_io, "persistent_run_i", action=OUTPUT)
      call print(persistent_run_i, file=persistent_io)
      call finalise(persistent_io)
   endif

   if (.not. persistent_already_started) then
      call print("writing POTCAR", PRINT_VERBOSE)
      call write_vasp_potcar(at, trim(run_dir), trim(potcar_files), nelect_of_elem, error=error)
      PASS_ERROR(error)
   endif

   ! total number of electrons for each element from POTCAR files
   nelect = sum(nelect_of_elem(at%Z))
   ! get Charge from xyz file
   r_charge = 0.0_dp
   if (get_value(at%params, 'Charge', r_charge)) then
      nelect = nelect - r_charge
   else if (get_value(at%params, 'Charge', i_charge)) then
      r_charge = i_charge
      nelect = nelect - i_charge
   endif
   if (r_charge .fne. 0.0_dp) then ! set total number of electrons
      call set_value(incar_dict, "NELECT", nelect)
   endif

   if (assign_pointer(at, "magmom", magmom_p)) then
      call set_value(incar_dict, "MAGMOM", magmom_p)
   endif

   if (.not. persistent_already_started) then
      ! write INCAR
      call print("writing INCAR", PRINT_VERBOSE)
      call initialise(io, trim(run_dir)//"/INCAR", action=OUTPUT, error=error)
      PASS_ERROR(error)
      call print(write_string(incar_dict, entry_sep=quip_new_line, quote_char=' '), file=io)
      call finalise(io)

      ! write KPOINTS if needed
      if (trim(kpoints_file) /= "_NONE_") then
         call print("writing KPOINTS", PRINT_VERBOSE)
         call system_command("cp "//trim(kpoints_file)//" "//trim(run_dir)//"/KPOINTS", status=stat)
         if (stat /= 0) then
            RAISE_ERROR("do_vasp_calc copying kpoints_file "//trim(kpoints_file)//" into "//trim(run_dir), error)
         endif
      endif
   endif

   if (.not. persistent_already_started .or. mod(persistent_run_i, persistent_restart_interval) == 0 .or. at%n /= persistent_n_atoms) then
      call print("writing POSCAR", PRINT_VERBOSE)
      call write_vasp_poscar(at, trim(run_dir), error=error)
      PASS_ERROR(error)
   endif
   call system_timer('vasp_filepot/prep')

   if (persistent) then ! restart VASP if needed, write to REFTRAJCAR if needed, initiate calc
      call system_timer('vasp_filepot/start_restart')
      if (mod(persistent_run_i, persistent_restart_interval) == 0 .or. at%n /= persistent_n_atoms) then ! start or restart VASP
         call print("starting or restarting VASP", PRINT_VERBOSE)
         if (persistent_already_started) then ! tell currently running vasp to quit
            call initialise(persistent_io, trim(run_dir)//"/REFTRAJCAR", action=OUTPUT)
            call print("0", file=persistent_io)
            call finalise(persistent_io)
            call initialise(persistent_io, trim(run_dir)//"/REFTRAJ_READY", action=OUTPUT)
            call print("-", file=persistent_io)
            call finalise(persistent_io)
            ! wait for signal that it's finished
            call system_command("echo 'before wait'; ps auxw | egrep 'vasp|mpi'")
            call wait_for_file_to_exist(trim(run_dir)//"/vasp.done", max_wait_time=600.0_dp, error=error)
            call system_command("echo 'after wait'; ls -l "//trim(run_dir))
            call system_command("echo 'after wait'; ps auxw | egrep 'vasp|mpi'")
            PASS_ERROR(error)
            ! clean up
            call unlink(trim(run_dir)//"/vasp.done")
            call system_command("cd "//trim(run_dir)//"; rm -f REFTRAJCAR REFTRAJ_READY OUTCAR CONTCAR OSZICAR")
            call system_command("cat vasp.stdout."//run_dir_i//" >> vasp.stdout")
         endif
         ! start vasp in background
         call          print("cd "//trim(run_dir)//"; bash -c "//'"'//"("//'\"'//trim(vasp_path)//'\"'//" > ../vasp.stdout."//run_dir_i//"; echo done > vasp.done )"//'"'//" 2>&1 &")
         call system_command("cd "//trim(run_dir)//"; bash -c "//'"'//"("//'\"'//trim(vasp_path)//'\"'//" > ../vasp.stdout."//run_dir_i//"; echo done > vasp.done )"//'"'//" 2>&1 &", status=stat)
         if (stat /= 0) then
            RAISE_ERROR("error running "//trim(vasp_path)//", status="//stat, error)
         endif
         persistent_already_started = .true.
      endif
      call system_timer('vasp_filepot/start_restart')
      call system_timer('vasp_filepot/write_reftrajcar')
      if (mod(persistent_run_i,persistent_restart_interval) /= 0 .and. at%n == persistent_n_atoms) then ! need to write current coords to REFTRAJCAR
         call print("writing current coords to REFTRAJCAR", PRINT_VERBOSE)
         call initialise(persistent_io, trim(run_dir)//"/REFTRAJCAR", action=OUTPUT)
         call print(at%N, file=persistent_io)
         call fix_lattice(at)
         call print(at%lattice(1,:), file=persistent_io)
         call print(at%lattice(2,:), file=persistent_io)
         call print(at%lattice(3,:), file=persistent_io)
         do i=1, at%N
            call print(matmul(at%g, at%pos(:,i)), file=persistent_io)
         end do
         call finalise(persistent_io)
         call print("touching REFTRAJ_READY", PRINT_VERBOSE)
         call initialise(persistent_io, trim(run_dir)//"/REFTRAJ_READY", action=OUTPUT)
         call print("-", file=persistent_io)
         call finalise(persistent_io)
      endif
      call system_timer('vasp_filepot/write_reftrajcar')
      call system_timer('vasp_filepot/wait_reftraj_step_done')
      call print("waiting for output REFTRAJ_STEP_DONE", PRINT_VERBOSE)
      ! wait for output to be ready
      call wait_for_file_to_exist(trim(run_dir)//"/REFTRAJ_STEP_DONE", persistent_max_wait_time, cycle_time=0.1_dp, error=error)
      PASS_ERROR(error)
      call system_timer('vasp_filepot/wait_reftraj_step_done')
   else ! not persistent, just run vasp
      call system_timer('vasp_filepot/run_vasp')
      call print("running vasp", PRINT_VERBOSE)
      call print("cd "//trim(run_dir)//"; bash -c "//'"'//trim(vasp_path)//'"'//" > ../vasp.stdout."//run_dir_i//" 2>&1")
      call system_command("cd "//trim(run_dir)//"; bash -c "//'"'//trim(vasp_path)//'"'//" > ../vasp.stdout."//run_dir_i//" 2>&1", status=stat)
      if (stat /= 0) then
         RAISE_ERROR("error running "//trim(vasp_path)//", status="//stat, error)
      endif
      call system_command("cat vasp.stdout."//run_dir_i//" >> vasp.stdout")
      call system_timer('vasp_filepot/run_vasp')
   endif

   ! look for errors in stdout and runtime in OUTCAR
   call system_command("fgrep -i 'error' vasp.stdout."//run_dir_i)
   call system_command("fgrep -i 'total cpu time used' "//trim(run_dir)//"/OUTCAR | tail -1; fgrep -c 'LOOP:' "//trim(run_dir)//"/OUTCAR | tail -1")

   call system_timer('vasp_filepot/read_output')
   if (persistent) then ! read REFTRAJ output
      call print("reading output REFTRAJ_OUTPUT", PRINT_VERBOSE)
      ! read output
      call read_vasp_output_persistent(at, run_dir, do_calc_energy, do_calc_force, do_calc_virial, error)
      PASS_ERROR(error)
      converged = .true. ! should really find a way of telling if converged for persistent connections
      ! clean up
      call unlink(trim(run_dir)//"/REFTRAJ_OUTPUT")
      call unlink(trim(run_dir)//"/REFTRAJ_STEP_DONE")
   else ! read regular VASP output (OUTCAR)
      call print("reading vasp output", PRINT_VERBOSE)
      call read_vasp_output(trim(run_dir), do_calc_energy, do_calc_force, do_calc_virial, converged, error=error)
      PASS_ERROR(error)
   endif
   call system_timer('vasp_filepot/read_output')

   call system_timer('vasp_filepot/end')
   if (.not. converged) then
      call print("WARNING: failed to find sign of convergence in OUTCAR", verbosity=PRINT_ALWAYS)
      if (.not. ignore_convergence) then
         RAISE_ERROR("failed to find sign of convergence in OUTCAR", error)
      endif
   endif

   inquire(file=trim(run_dir)//"/WAVECAR", exist=have_wavecar)
   if (have_wavecar) then
      call print("copying run_dir/WAVECAR to WAVECAR.run_suffix", PRINT_VERBOSE)
      call system_command("cp "//trim(run_dir)//"/WAVECAR WAVECAR."//trim(run_suffix), status=stat)
      if (stat /= 0) then
         RAISE_ERROR("do_vasp_calc copying file "//trim(run_dir)//"WAVECAR into ./"//trim(run_dir)//"WAVECAR."//trim(run_suffix), error)
      endif
   endif

   call print("doing unsorting", PRINT_VERBOSE)
   call atoms_sort(at, 'sort_index', error=error)
   PASS_ERROR(error)

   if (clean_up_files .and. .not. persistent) then
      call print("cleaning up", PRINT_VERBOSE)
      call system_command("rm -rf vasp_run_"//(1+mod(run_dir_i,clean_up_keep_n+1)))
   endif

   call system_timer('vasp_filepot/end')

   call verbosity_pop()

end subroutine do_vasp_calc

subroutine write_vasp_potcar(at, run_dir, potcar_files, nelect_of_elem, error)
   type(Atoms), intent(in) :: at
   character(len=*), intent(in) :: run_dir, potcar_files
   real(dp), intent(out) :: nelect_of_elem(:)
   integer, intent(out), optional :: error

   integer :: i, Z_i, potcar_files_n_fields
   character(len=STRING_LENGTH) :: potcar_files_fields(128), potcar_files_a(128)
   integer, allocatable :: uniq_Z(:)

   integer :: stat, t_Z
   real(dp) :: t_nelect
   character(len=STRING_LENGTH) :: line
   type(inoutput) :: io

   INIT_ERROR(error)

   call split_string_simple(trim(potcar_files), potcar_files_fields, potcar_files_n_fields, " ", error)
   PASS_ERROR(error)

   if (potcar_files_n_fields <= 0 .or. mod(potcar_files_n_fields,2) /= 0) then
      RAISE_ERROR("write_vasp_potcar didn't find even number of fields (Z1 potcar1 [Z2 potcar2 ...]) in potcar_files='"//trim(potcar_files)//"'", error)
   endif

   call uniq(at%Z, uniq_Z)

   potcar_files_a=""
   do i=1, potcar_files_n_fields, 2
      read(unit=potcar_files_fields(i),fmt=*) Z_i
      potcar_files_a(Z_i) = trim(potcar_files_fields(i+1))
   end do

   call system_command("rm -f nelect_summary; touch nelect_summary")
   call system_command("rm -f "//trim(run_dir)//"/POTCAR")
   do Z_i=1, size(uniq_Z)
      if (uniq_Z(Z_i) == 0) exit
      if (len_trim(potcar_files_a(uniq_Z(Z_i))) == 0) then
         RAISE_ERROR("write_vasp_potcar could not find a potcar for Z="//uniq_Z(Z_i), error)
      endif
      call system_command("cat "//trim(potcar_files_a(uniq_Z(Z_i)))//" >> "//trim(run_dir)//"/POTCAR")
      call system_command("(echo -n '"//uniq_Z(Z_i)//" '; cat "//trim(potcar_files_a(uniq_Z(Z_i)))//" | fgrep ZVAL | sed -e 's/.*ZVAL[ ]*=[ ]*//' -e 's/ .*//') >> nelect_summary")
   end do

   ! read back nelect_summary
   call initialise(io, "nelect_summary")
   do Z_i=1, size(uniq_Z)
      line = read_line(io, stat)
      read (unit=line, fmt=*) t_Z, t_nelect
      if (t_Z <= 0 .or. t_Z > size(nelect_of_elem)) then
         RAISE_ERROR("write_vasp_potcar got value out of range reading nelect_summary value t_Z="//t_Z, error)
      endif
      if (t_Z /= uniq_Z(Z_i)) then
         RAISE_ERROR("write_vasp_potcar got mismatch reading nelect_summary value t_Z="//t_Z//" uniq_Z("//Z_i//")="//uniq_Z(Z_i), error)
      endif
      nelect_of_elem(t_Z) = t_nelect
   end do
   call finalise(io)

end subroutine write_vasp_potcar

subroutine write_vasp_poscar(at, run_dir, error)
   type(Atoms), intent(inout) :: at
   character(len=*), intent(in) :: run_dir
   integer, intent(out), optional :: error

   integer, allocatable :: uniq_Z(:), n_Z(:)
   integer :: i
   logical :: swapped_a1_a2

   type(inoutput) :: io

   INIT_ERROR(error)

   call fix_lattice(at, swapped_a1_a2)

   call initialise(io, trim(run_dir)//"/POSCAR",action=OUTPUT)
   if (swapped_a1_a2) then
      call print("QUIP VASP run "//trim(run_dir)//" swapped a1 and a2", file=io)
   else
      call print("QUIP VASP run "//trim(run_dir), file=io)
   endif
   call print("1.00 ! scale factor", file=io)

   call print(at%lattice(:,1), file=io)
   call print(at%lattice(:,2), file=io)
   call print(at%lattice(:,3), file=io)

   call uniq(at%Z, uniq_Z)
   allocate(n_Z(size(uniq_Z)))
   do i=1, size(uniq_Z)
      n_Z(i) = count(at%Z == uniq_Z(i))
   end do
   call print(""//n_Z//" ! Zs = "//uniq_Z, file=io)
   call print("Selective Dynamics", file=io)
   call print("Cartesian", file=io)
   do i=1, at%N
      call print(at%pos(:,i)//" T T T", file=io)
   end do
   call finalise(io)

end subroutine write_vasp_poscar

subroutine fix_lattice(at, swapped_a1_a2)
   type(Atoms), intent(inout) :: at
   logical, intent(out), optional :: swapped_a1_a2

   real(dp) :: vol, lattice(3,3)

   if (present(swapped_a1_a2)) swapped_a1_a2 = .false.

   vol = scalar_triple_product(at%lattice(:,1), at%lattice(:,2), at%lattice(:,3))
   if (vol < 0.0_dp) then
      if (present(swapped_a1_a2)) swapped_a1_a2 = .true.
      lattice(:,1) = at%lattice(:,2)
      lattice(:,2) = at%lattice(:,1)
      lattice(:,3) = at%lattice(:,3)
      call set_lattice(at, lattice, .false.)
   endif

end subroutine fix_lattice

subroutine read_vasp_output(run_dir, do_calc_energy, do_calc_force, do_calc_virial, converged, error)
   character(len=*), intent(in) :: run_dir
   logical, intent(in) :: do_calc_energy, do_calc_force, do_calc_virial
   logical, intent(out) :: converged
   integer, intent(out), optional :: ERROR

   type(inoutput) :: outcar_io
   integer :: stat
   character(len=STRING_LENGTH) :: line
   character(len=STRING_LENGTH)  :: fields(100), t_s
   integer :: n_fields, line_i
   logical in_virial

   real(dp) :: energy, virial(3,3), t_pos(3)
   real(dp), pointer :: force_p(:,:)
   integer :: force_start_line_i, virial_start_line_i, read_err

   INIT_ERROR(error)

   if (do_calc_force) then
      if (.not. assign_pointer(at, "force", force_p)) then
         call add_property(at, "force", 0.0_dp, n_cols=3, ptr2=force_p, error=error)
         PASS_ERROR(error)
      endif
   endif

   converged = .false.
   force_start_line_i = -1
   in_virial = .false.
   call initialise(outcar_io, trim(run_dir)//"/OUTCAR")
   stat = 0
   line_i = 1
   do while (stat == 0)
      line = read_line(outcar_io, stat)
      if (do_calc_energy) then
         if (index(trim(line),"free  energy") > 0) then ! found energy
            call split_string_simple(trim(line), fields, n_fields, " ", error=error)
            PASS_ERROR(error)
            if (n_fields /= 6) then
               RAISE_ERROR("read_vasp_output confused by energy line# "//line_i//" in OUTCAR line '"//trim(line)//"'", error)
            endif
            if (fields(6) /= "eV") then
               RAISE_ERROR("read_vasp_output confused by units energy line# "//line_i//" in OUTCAR line '"//trim(line)//"'", error)
            endif
            read(unit=fields(5),fmt=*) energy
            call set_param_value(at, 'energy', energy)
         endif
      endif
      if (do_calc_force) then
         if (force_start_line_i > 0 .and. line_i >= force_start_line_i .and. line_i <= force_start_line_i+at%N-1) then
            read(unit=line,fmt=*) t_pos, force_p(:,line_i-force_start_line_i+1)
         else
            if (index(trim(line),"TOTAL-FORCE") > 0) then ! found force
               call split_string_simple(trim(line), fields, n_fields, " ", error=error)
               if (fields(3) /= "(eV/Angst)") then
                  RAISE_ERROR("read_vasp_output confused by units force header line# "//line_i//" in OUTCAR line '"//trim(line)//"'", error)
               endif
               force_start_line_i = line_i+2
            endif
         endif
      endif
      if (do_calc_virial) then
         if (in_virial .and. index(trim(line),"Total") > 0) then
            read(unit=line,fmt=*,iostat=read_err) t_s, virial(1,1), virial(2,2), virial(3,3), virial(1,2), virial(2,3), virial(1,3)
            if (read_err /= 0) then
               RAISE_ERROR ("read_vasp_output failed to read virial from line '"//trim(line)//"'", error)
            endif
            virial(2,1) = virial(1,2)
            virial(3,2) = virial(2,3)
            virial(3,1) = virial(1,3)
            call set_param_value(at, 'virial', virial)
            in_virial = .false.
         else
            if (index(trim(line),"FORCE on cell") > 0) then ! found start of virial section
               call split_string_simple(trim(line), fields, n_fields, " ", error=error)
               if (fields(9) /= "(eV):" .and. fields(9) /= "(eV/reduce") then
                  RAISE_ERROR("read_vasp_output confused by units virial header line# "//line_i//" in OUTCAR line '"//trim(line)//"'", error)
               endif
               in_virial = .true.
            endif
         endif
      endif
      if (.not. converged) then
         if (index(trim(line),"aborting loop because EDIFF is reached") > 0) then ! found convergence of SCF
            converged = .true.
         endif
      endif
      ! if (stat == 0) call print("GOT OUTPUT "//trim(line))
      line_i = line_i + 1
   end do
   call finalise(outcar_io)

end subroutine read_vasp_output

subroutine read_vasp_output_persistent(at, run_dir, do_calc_energy, do_calc_force, do_calc_virial, error)
   type(atoms), intent(inout) :: at
   character(len=*), intent(in) :: run_dir
   logical, intent(in) :: do_calc_energy, do_calc_force, do_calc_virial
   integer, optional, intent(out) :: error

   real(dp) :: energy, virial(3,3), force_t(3)
   real(dp), pointer :: force_p(:,:)
   type(inoutput) :: reftraj_output_io
   character(len=STRING_LENGTH) :: line
   integer :: ni

   INIT_ERROR(error)

   if (do_calc_force) then
      if (.not. assign_pointer(at, "force", force_p)) then
         call add_property(at, "force", 0.0_dp, n_cols=3, ptr2=force_p, error=error)
         PASS_ERROR(error)
      endif
   endif

   call initialise(reftraj_output_io, trim(run_dir)//"/REFTRAJ_OUTPUT", action=INPUT)

   line = read_line(reftraj_output_io); read(unit=line, fmt=*) ni
   if (ni /= at%N) then
      RAISE_ERROR("read_vasp_output_persistent got NI == "//ni//" != at%N == "//at%N, error)
   endif

   ! read energy
   line = read_line(reftraj_output_io); read(unit=line, fmt=*) energy
   if (do_calc_energy) call set_param_value(at, 'energy', energy)
   ! read force
   do i=1, at%N
      line = read_line(reftraj_output_io); read(unit=line, fmt=*) force_t
      if (do_calc_force) force_p(:,i) = force_t(:)
   end do
   ! read virial
   if (do_calc_virial) then
      line = read_line(reftraj_output_io); read(unit=line, fmt=*) virial(1,1), virial(2,2), virial(3,3), virial(1,2), virial(2,3), virial(3,1)
      virial(2,1) = virial(1,2)
      virial(3,2) = virial(2,3)
      virial(1,3) = virial(3,1)
      call set_param_value(at, "virial", virial)
   endif

   call finalise(reftraj_output_io)
end subroutine read_vasp_output_persistent

subroutine read_vasp_incar_dict(incar_dict, incar_template_file, error)
   type(Dictionary), intent(inout) :: incar_dict
   character(len=*), intent(in) :: incar_template_file
   integer, intent(out), optional :: error

   type(Inoutput) :: incar_io
   character(len=STRING_LENGTH), allocatable :: incar_a(:), incar_line_fields(:), incar_field_fields(:)
   integer :: incar_n_lines, incar_line_n_fields, incar_field_n_fields
   integer :: i, j
   integer :: comment_pos

   INIT_ERROR(error)

   call initialise(incar_io, trim(incar_template_file), INPUT)
   call read_file(incar_io, incar_a, incar_n_lines)
   call finalise(incar_io)

   call initialise(incar_dict)
   allocate(incar_line_fields(100))
   allocate(incar_field_fields(100))
   ! loop over lines
   do i=1, incar_n_lines
      comment_pos = index(trim(incar_a(i)), '#')
      if (comment_pos > 1) then
         incar_a(i) = incar_a(i)(1:comment_pos-1)
      elseif (comment_pos == 1) then
         incar_a(i) = ""
      endif
      comment_pos = index(trim(incar_a(i)), '!')
      if (comment_pos > 1) then
         incar_a(i) = incar_a(i)(1:comment_pos-1)
      elseif (comment_pos == 1) then
         incar_a(i) = ""
      endif
      call split_string(trim(incar_a(i)), ";", "''"//'""', incar_line_fields, incar_line_n_fields, .true.)
      ! loop over ';' separated fields
      do j=1, incar_line_n_fields
         call split_string(incar_line_fields(j), "=", "''"//'""', incar_field_fields, incar_field_n_fields, .true.)
         if (incar_field_n_fields > 2) then
            RAISE_ERROR("read_vasp_incar_dict got more than one '=' in field "//trim(incar_line_fields(j)), error)
         endif
         if (.not. has_key(incar_dict, trim(adjustl(incar_field_fields(1))))) then
            ! VASP always uses first instance of a keyword, do same here
            call set_value(incar_dict, trim(adjustl(incar_field_fields(1))), trim(adjustl(incar_field_fields(2))))
         endif
      end do
   end do

   deallocate(incar_line_fields)
   deallocate(incar_field_fields)

end subroutine read_vasp_incar_dict

end program vasp_driver

max_stars_repo_name: libAtoms/QUIP

max_issues_repo_name: libAtoms/QUIP

max_forks_repo_name: libAtoms/QUIP


text:
*DECK XPQNU
      SUBROUTINE XPQNU (NU1, NU2, MU, THETA, ID, PQA, IPQA, IERROR)
C***BEGIN PROLOGUE  XPQNU
C***SUBSIDIARY
C***PURPOSE  To compute the values of Legendre functions for XLEGF.
C            This subroutine calculates initial values of P or Q using
C            power series, then performs forward nu-wise recurrence to
C            obtain P(-MU,NU,X), Q(0,NU,X), or Q(1,NU,X). The nu-wise
C            recurrence is stable for P for all mu and for Q for mu=0,1.
C***LIBRARY   SLATEC
C***CATEGORY  C3A2, C9
C***TYPE      SINGLE PRECISION (XPQNU-S, DXPQNU-D)
C***KEYWORDS  LEGENDRE FUNCTIONS
C***AUTHOR  Smith, John M., (NBS and George Mason University)
C***ROUTINES CALLED  XADD, XADJ, XPSI
C***COMMON BLOCKS    XBLK1
C***REVISION HISTORY  (YYMMDD)
C   820728  DATE WRITTEN
C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
C   901019  Revisions to prologue.  (DWL and WRB)
C   901106  Changed all specific intrinsics to generic.  (WRB)
C           Corrected order of sections in prologue and added TYPE
C           section.  (WRB)
C   920127  Revised PURPOSE section of prologue.  (DWL)
C***END PROLOGUE  XPQNU
      REAL A,NU,NU1,NU2,PQ,PQA,XPSI,R,THETA,W,X,X1,X2,XS,
     1 Y,Z
      REAL DI,DMU,PQ1,PQ2,FACTMU,FLOK
      DIMENSION PQA(*),IPQA(*)
      COMMON /XBLK1/ NBITSF
      SAVE /XBLK1/
C
C        J0, IPSIK, AND IPSIX ARE INITIALIZED IN THIS SUBROUTINE.
C        J0 IS THE NUMBER OF TERMS USED IN SERIES EXPANSION
C        IN SUBROUTINE XPQNU.
C        IPSIK, IPSIX ARE VALUES OF K AND X RESPECTIVELY
C        USED IN THE CALCULATION OF THE XPSI FUNCTION.
C
C***FIRST EXECUTABLE STATEMENT  XPQNU
      IERROR=0
      J0=NBITSF
      IPSIK=1+(NBITSF/10)
      IPSIX=5*IPSIK
      IPQ=0
C        FIND NU IN INTERVAL [-.5,.5) IF ID=2  ( CALCULATION OF Q )
      NU=MOD(NU1,1.)
      IF(NU.GE..5) NU=NU-1.
C        FIND NU IN INTERVAL (-1.5,-.5] IF ID=1,3, OR 4  ( CALC. OF P )
      IF(ID.NE.2.AND.NU.GT.-.5) NU=NU-1.
C        CALCULATE MU FACTORIAL
      K=MU
      DMU=MU
      IF(MU.LE.0) GO TO 60
      FACTMU=1.
      IF=0
      DO 50 I=1,K
      FACTMU=FACTMU*I
   50 CALL XADJ(FACTMU,IF,IERROR)
      IF (IERROR.NE.0) RETURN
   60 IF(K.EQ.0) FACTMU=1.
      IF(K.EQ.0) IF=0
C
C        X=COS(THETA)
C        Y=SIN(THETA/2)**2=(1-X)/2=.5-.5*X
C        R=TAN(THETA/2)=SQRT((1-X)/(1+X)
C
      X=COS(THETA)
      Y=SIN(THETA/2.)**2
      R=TAN(THETA/2.)
C
C        USE ASCENDING SERIES TO CALCULATE TWO VALUES OF P OR Q
C        FOR USE AS STARTING VALUES IN RECURRENCE RELATION.
C
      PQ2=0.0
      DO 100 J=1,2
      IPQ1=0
      IF(ID.EQ.2) GO TO 80
C
C        SERIES FOR P ( ID = 1, 3, OR 4 )
C        P(-MU,NU,X)=1./FACTORIAL(MU)*SQRT(((1.-X)/(1.+X))**MU)
C                *SUM(FROM 0 TO J0-1)A(J)*(.5-.5*X)**J
C
      IPQ=0
      PQ=1.
      A=1.
      IA=0
      DO 65 I=2,J0
      DI=I
      A=A*Y*(DI-2.-NU)*(DI-1.+NU)/((DI-1.+DMU)*(DI-1.))
      CALL XADJ(A,IA,IERROR)
      IF (IERROR.NE.0) RETURN
      IF(A.EQ.0.) GO TO 66
      CALL XADD(PQ,IPQ,A,IA,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
   65 CONTINUE
   66 CONTINUE
      IF(MU.LE.0) GO TO 90
      X2=R
      X1=PQ
      K=MU
      DO 77 I=1,K
      X1=X1*X2
   77 CALL XADJ(X1,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      PQ=X1/FACTMU
      IPQ=IPQ-IF
      CALL XADJ(PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      GO TO 90
C
C        Z=-LN(R)=.5*LN((1+X)/(1-X))
C
   80 Z=-LOG(R)
      W=XPSI(NU+1.,IPSIK,IPSIX)
      XS=1./SIN(THETA)
C
C        SERIES SUMMATION FOR Q ( ID = 2 )
C        Q(0,NU,X)=SUM(FROM 0 TO J0-1)((.5*LN((1+X)/(1-X))
C    +XPSI(J+1,IPSIK,IPSIX)-XPSI(NU+1,IPSIK,IPSIX)))*A(J)*(.5-.5*X)**J
C
C        Q(1,NU,X)=-SQRT(1./(1.-X**2))+SQRT((1-X)/(1+X))
C             *SUM(FROM 0 T0 J0-1)(-NU*(NU+1)/2*LN((1+X)/(1-X))
C                 +(J-NU)*(J+NU+1)/(2*(J+1))+NU*(NU+1)*
C     (XPSI(NU+1,IPSIK,IPSIX)-XPSI(J+1,IPSIK,IPSIX))*A(J)*(.5-.5*X)**J
C
C        NOTE, IN THIS LOOP K=J+1
C
      PQ=0.
      IPQ=0
      IA=0
      A=1.
      DO 85 K=1,J0
      FLOK=K
      IF(K.EQ.1) GO TO 81
      A=A*Y*(FLOK-2.-NU)*(FLOK-1.+NU)/((FLOK-1.+DMU)*(FLOK-1.))
      CALL XADJ(A,IA,IERROR)
      IF (IERROR.NE.0) RETURN
   81 CONTINUE
      IF(MU.GE.1) GO TO 83
      X1=(XPSI(FLOK,IPSIK,IPSIX)-W+Z)*A
      IX1=IA
      CALL XADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      GO TO 85
   83 X1=(NU*(NU+1.)*(Z-W+XPSI(FLOK,IPSIK,IPSIX))+(NU-FLOK+1.)
     1  *(NU+FLOK)/(2.*FLOK))*A
      IX1=IA
      CALL XADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
   85 CONTINUE
      IF(MU.GE.1) PQ=-R*PQ
      IXS=0
      IF(MU.GE.1) CALL XADD(PQ,IPQ,-XS,IXS,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      IF(J.EQ.2) MU=-MU
      IF(J.EQ.2) DMU=-DMU
   90 IF(J.EQ.1) PQ2=PQ
      IF(J.EQ.1) IPQ2=IPQ
      NU=NU+1.
  100 CONTINUE
      K=0
      IF(NU-1.5.LT.NU1) GO TO 120
      K=K+1
      PQA(K)=PQ2
      IPQA(K)=IPQ2
      IF(NU.GT.NU2+.5) RETURN
  120 PQ1=PQ
      IPQ1=IPQ
      IF(NU.LT.NU1+.5) GO TO 130
      K=K+1
      PQA(K)=PQ
      IPQA(K)=IPQ
      IF(NU.GT.NU2+.5) RETURN
C
C        FORWARD NU-WISE RECURRENCE FOR F(MU,NU,X) FOR FIXED MU
C        USING
C        (NU+MU+1)*F(MU,NU,X)=(2.*NU+1)*F(MU,NU,X)-(NU-MU)*F(MU,NU-1,X)
C        WHERE F(MU,NU,X) MAY BE P(-MU,NU,X) OR IF MU IS REPLACED
C        BY -MU THEN F(MU,NU,X) MAY BE Q(MU,NU,X).
C        NOTE, IN THIS LOOP, NU=NU+1
C
  130 X1=(2.*NU-1.)/(NU+DMU)*X*PQ1
      X2=(NU-1.-DMU)/(NU+DMU)*PQ2
      CALL XADD(X1,IPQ1,-X2,IPQ2,PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      CALL XADJ(PQ,IPQ,IERROR)
      IF (IERROR.NE.0) RETURN
      NU=NU+1.
      PQ2=PQ1
      IPQ2=IPQ1
      GO TO 120
C
      END


text:
      PROGRAM DSEX04
C
      PARAMETER (NUM=16, NX=21, NY=21)
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
      REAL XI(NUM), YI(NUM), ZI(NUM)
      REAL XO(NX), YO(NY), OUTPUT(NX,NY)
      REAL XINC, YINC
C
C  Input data points and values.
C
      DATA XI/0.00, 1.00, 0.00, 1.00, 0.30, 0.30, 0.30, 0.69,
     +        0.71, 0.71, 0.69, 0.70, 0.70, 0.70, 0.69, 0.71/
      DATA YI/0.00, 0.00, 1.00, 1.00, 0.70, 0.30, 0.70, 0.69,
     +        0.69, 0.71, 0.71, 0.70, 0.69, 0.71, 0.70, 0.70/
      DATA ZI/0.00, 0.00, 0.00, 0.50, 0.50, 0.50, 0.50, 1.00,
     +        1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00/
      DATA RHO, THETA, PHI/3.0, -45., 55./
C
C  Specify the output grid.
C
      XINC = 1./ REAL(NX-1)
      YINC = 1./ REAL(NY-1)
      DO 30 I=1,NX
        XO(I) = REAL((I-1)*XINC)
        DO 40 J=1,NY
          YO(J) = REAL((J-1)*YINC)
   40   CONTINUE
   30 CONTINUE
C
C  Set shadowing flag.
C
      CALL DSSETI('SHD',1)
      CALL DSGRID2S(NUM, XI, YI, ZI, NX, NY, XO, YO, OUTPUT, IER)
      IF (IER .NE. 0) THEN
        WRITE(6,520) IER
  520   FORMAT(' Error ',I3,' returned from DSGRID2S')
        STOP
      ENDIF
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
      CALL TDEZ2D(NX, NY, XO, YO, OUTPUT, RHO, THETA, PHI, 6)
      CALL FRAME()
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
C
      STOP
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
subroutine init_random_seed()
  use iso_fortran_env, only: int64
  implicit none
  integer, allocatable :: seed(:)
  integer :: i, n, un, istat, dt(8), pid
  integer(int64) :: t

  call random_seed(size = n)
  allocate(seed(n))
  ! First try if the OS provides a random number generator
  open(newunit=un, file="/dev/urandom", access="stream", &
       form="unformatted", action="read", status="old", iostat=istat)
  if (istat == 0) then
     read(un) seed
     close(un)
  else
     ! Fallback to XOR:ing the current time and pid. The PID is
     ! useful in case one launches multiple instances of the same
     ! program in parallel.
     call system_clock(t)
     if (t == 0) then
        call date_and_time(values=dt)
        t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
             + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
             + dt(3) * 24_int64 * 60 * 60 * 1000 &
             + dt(5) * 60 * 60 * 1000 &
             + dt(6) * 60 * 1000 + dt(7) * 1000 &
             + dt(8)
     end if
     pid = getpid()
     t = ieor(t, int(pid, kind(t)))
     do i = 1, n
        seed(i) = lcg(t)
     end do
  end if
  call random_seed(put=seed)
contains
  ! This simple PRNG might not be good enough for real work, but is
  ! sufficient for seeding a better PRNG.
  function lcg(s)
    integer :: lcg
    integer(int64) :: s
    if (s == 0) then
       s = 104729
    else
       s = mod(s, 4294967296_int64)
    end if
    s = mod(s * 279470273_int64, 4294967291_int64)
    lcg = int(mod(s, int(huge(0), int64)), kind(0))
  end function lcg
end subroutine init_random_seed

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible