Merge pull request #5 from hokru/ciopt_patch

[dev branch] Enabled general external interface (GEI) for CIopt

docs/manual.rst

@@ -14,11 +14,12 @@ and allows code re-use and sharing as long as credits are given.
 Much of the development has been done in my free time and I learned a lot with this project.
 If you find bugs or have suggestions please write to me or open an issue on github.
 
+If you publish work using this code plesae be so kind to notify about the reference. It will help me secure future funding.
 
 How to cite
 ***********
 
-For proper technical documentation cite the git version, github repository and authors, e.g.:
+There is not yet any scientific publication for xopt. However, for proper technical documentation please cite the git version, github repository and authors. The version can be found in the output, e.g.:
 
 ``H.Kruse, https://github.com/hokru/xopt, version v2.0b-1-g7434``
 

docs/prog_interface.rst

@@ -75,7 +75,7 @@ following format in a.u.::
 
 
 
-The default name of the script is ``mygrad.sh`` and the default gradient file is ``xopt.grad``
+The default name of the script is ``mygrad.sh`` and the default gradient file is ``xopt.grad``. The script name can be changed either in the .xoptrc options or with ``-gei_command <script name>``.
 
 
 
@@ -151,7 +151,7 @@ We can do ``Gaussian``, ``Turbomole(dscf,ricc2)``, ``Orca`` and (modified) ``Amb
 
 You need to make 2 directories named ``stateJ.xopt` and ``stateI.xopt``.
 Prepare the input for each state inside the directories.
-It should work for Turbomole, ORCA and G09 if you follow the general preparation guidelines above.
+It should work for Turbomole(dscf/ricc2), ORCA, G09, AMBER and the GEI custom interface. For the GEI each subdirectory must contain it's own script that writes the ``xopt.grad`` file.
 
 We assume state I < state J, e.g. J=I+1.
 

extralibs/make.in

@@ -1,3 +1,3 @@
-#FC = gfortran -ffree-line-length-none -static 
+#FC = gfortran -ffree-line-length-none #-static 
 FC = ifort -static 
 FFLAGS= -O 

src/ciopt.f90

@@ -40,95 +40,102 @@ subroutine cioptgrad
     endif
 else
 
-
-! directory definitions. Fixed for now
-dir1='/stateI.xopt'
-dir2='/stateJ.xopt'
-!call getcwd(pwd)
-call get_environment_variable('PWD',pwd)
-
-! change to dir1
-exepath=adjustl(trim(pwd))//trim(adjustl(dir1))
-call chdir(adjustl(trim(exepath)))
-
-!----------------------------------------
-! for state I
-if (tm) then
-    call system('actual -r &> /dev/null')
-    aa='dscf > '//xjob
-    bb='ricc2 >> '//xjob
-    call system(aa)
-    call system(bb)
-    call tmgrad(nat,gi,ei)
-elseif(orca) then
-    aa='orca30 orca.in > '//xjob
-    call system(trim(aa))
-    call rdorcagrad(nat,gi,ei   ,'orca.engrad')
-    call system("cp xopt.job xopt.last")
-elseif(gaus) then
-    call IOgaussian('g.in')
-    aa='run-g09 g.in xopt.job'
-    call system(trim(aa))
-    call g09grad(nat,gi,ei)
-    call system("cp xopt.job xopt.last")
-elseif(amber) then
-    aa='rm -f forcedump.dat force.dat > /dev/null'
-    call system (trim(aa))
-    call wamber(nat3,xyz,'amber.rst')
-    call system(trim(command_amber))
-    call ambergrad(nat,gi,ei,apbs)
-    call cpfile('xopt.job','xopt.last','.')
-!    aa='sander -O -i '//trim(ambin)// &
-!             ' -o '//trim(xjob)// &
-!             ' -c amber.rst'// &
-!             ' -p '//trim(ambtop)
-!    call system(trim(aa))
-!    call ambergrad(nat,grad,energy,apbs)
-endif
-!----------------------------------------
-
-! change to dir2
-exepath=adjustl(trim(pwd))//trim(adjustl(dir2))
-call chdir(trim(exepath))
-
-!----------------------------------------
-! for state J
-if(tm) then
-    call system('actual -r &> /dev/null')
-    aa='dscf > '//xjob
-    bb='ricc2 >> '//xjob
-    call system(aa)
-    call system(bb)
-    call tmgrad(nat,gj,ej)
-elseif(orca) then
-    aa='orca30 orca.in > '//xjob
-    call system(trim(aa))
-    call rdorcagrad(nat,gj,ej,'orca.engrad')
-    call system("cp xopt.job xopt.last")
-elseif(gaus) then
-call IOgaussian('g.in')
-    aa='run-g09 g.in xopt.job'
-    call system(trim(aa))
-    call g09grad(nat,gj,ej)
-    call system("cp xopt.job xopt.last")
-elseif(amber) then
-    aa='rm -f forcedump.dat force.dat > /dev/null'
-    call system (trim(aa))
-    call wamber(nat3,xyz,'amber.rst')
-    call system(trim(command_amber))
-    call ambergrad(nat,gj,ej,apbs)
-    call cpfile('xopt.job','xopt.last','.')
-!    aa='sander -O -i '//trim(ambin)// &
-!             ' -o '//trim(xjob)// &
-!             ' -c amber.rst'// &
-!             ' -p '//trim(ambtop)
-!    call system(trim(aa))
-!    call ambergrad(nat,grad,energy,apbs)
-endif
-!----------------------------------------
-
-! return to working dir
-call chdir(adjustl(trim(pwd)))
+  ! directory definitions. Fixed for now
+  dir1='/stateI.xopt'
+  dir2='/stateJ.xopt'
+  !call getcwd(pwd)
+  call get_environment_variable('PWD',pwd)
+
+  ! change to dir1
+  exepath=adjustl(trim(pwd))//trim(adjustl(dir1))
+  call chdir(adjustl(trim(exepath)))
+
+  !----------------------------------------
+  ! It is an unfortunate design to have to repeat getgrad.
+  ! Long term this should be done better
+  ! for state I
+  if(gei) then
+    call system(trim(command_gei))
+    call egradfile(nat,gi,ei,trim(geigrad))
+  elseif (tm) then
+      call system('actual -r &> /dev/null')
+      aa='dscf > '//xjob
+      bb='ricc2 >> '//xjob
+      call system(aa)
+      call system(bb)
+      call tmgrad(nat,gi,ei)
+  elseif(orca) then
+      call system(trim(command_orca))
+      call rdorcagrad(nat,gi,ei   ,'orca.engrad')
+      call system("cp xopt.job xopt.last")
+  elseif(gaus) then
+      call IOgaussian('g.in')
+      call system(trim(command_gaus))
+      call g09grad(nat,gi,ei)
+      call system("cp xopt.job xopt.last")
+  elseif(amber) then
+      aa='rm -f forcedump.dat force.dat > /dev/null'
+      call system (trim(aa))
+      call wamber(nat3,xyz,'amber.rst')
+      call system(trim(command_amber))
+      call ambergrad(nat,gi,ei,apbs)
+      call cpfile('xopt.job','xopt.last','.')
+  !    aa='sander -O -i '//trim(ambin)// &
+  !             ' -o '//trim(xjob)// &
+  !             ' -c amber.rst'// &
+  !             ' -p '//trim(ambtop)
+  !    call system(trim(aa))
+  !    call ambergrad(nat,grad,energy,apbs)
+  endif
+  ! call getd3gcp(nat,ei,gi)
+  !----------------------------------------
+
+  ! change to dir2
+  exepath=adjustl(trim(pwd))//trim(adjustl(dir2))
+  call chdir(trim(exepath))
+
+  !----------------------------------------
+  ! for state J
+  if(gei) then
+    call system(trim(command_gei))
+    call egradfile(nat,gj,ej,trim(geigrad))
+  elseif (tm) then
+      call system('actual -r &> /dev/null')
+      aa='dscf > '//xjob
+      bb='ricc2 >> '//xjob
+      call system(aa)
+      call system(bb)
+      call tmgrad(nat,gj,ej)
+  elseif(orca) then
+      aa='orca30 orca.in > '//xjob
+      call system(trim(aa))
+      call rdorcagrad(nat,gj,ej,'orca.engrad')
+      call system("cp xopt.job xopt.last")
+  elseif(gaus) then
+      call IOgaussian('g.in')
+      aa='run-g09 g.in xopt.job'
+      call system(trim(aa))
+      call g09grad(nat,gj,ej)
+      call system("cp xopt.job xopt.last")
+  elseif(amber) then
+      aa='rm -f forcedump.dat force.dat > /dev/null'
+      call system (trim(aa))
+      call wamber(nat3,xyz,'amber.rst')
+      call system(trim(command_amber))
+      call ambergrad(nat,gj,ej,apbs)
+      call cpfile('xopt.job','xopt.last','.')
+  !    aa='sander -O -i '//trim(ambin)// &
+  !             ' -o '//trim(xjob)// &
+  !             ' -c amber.rst'// &
+  !             ' -p '//trim(ambtop)
+  !    call system(trim(aa))
+  !    call ambergrad(nat,grad,energy,apbs)
+  endif
+  ! call getd3gcp(nat,ej,gj)
+  !----------------------------------------
+
+  ! return to working dir
+  call chdir(adjustl(trim(pwd)))
 endif
 
 flipped=.false.

src/eval_opt.f90

@@ -111,7 +111,7 @@ subroutine loadrc()
 subroutine eval_options (infile)
 use logic
 use popt
-use progs, only: usrscr,prog_flags
+use progs, only: usrscr,prog_flags, scall_gei
 use MDdat
 use fiso, only: r8, stdout,stdout_default
 implicit none
@@ -177,6 +177,9 @@ subroutine eval_options (infile)
     gei=.true.
 !    call_gei=trim(narg) ??
   endif
+  if(index(ftmp,'-gei_command ').ne.0) then
+    scall_gei=trim(adjustl(arg(i+1)))
+  endif
   if(index(ftmp,'-ppot ').ne.0) then
    ppot=.true.
    tm=.true.
@@ -441,6 +444,7 @@ subroutine help
 !write(io,'(a)')''-noprint
 write(io,'(a)')''
 write(io,'(a)')'PROGRAM INTERFACES'
+ write(io,'(a)')''
  write(io,'(a)')'-tm                Turbomole(Xscf/Xgrad)'
  write(io,'(a)')'-tmcc              Turbomole(Xscf/ricc2)'
  write(io,'(a)')'-grad              stop after gradient(includes restraints!)'
@@ -466,6 +470,10 @@ subroutine help
  write(io,'(a)')''
  write(io,'(a)')'-numgrad           arbitrary numerical gradient (provide mygrad.sh)'
  write(io,'(a)')''
+ write(io,'(a)')'CUSTOM INTERFACE'
+ write(io,'(a)')'-gei               use general external interface'
+ write(io,'(a)')'-gei_command       name of the script to call'
+ write(io,'(a)')''
  write(io,'(a)')'HESSIANS'
  write(io,'(a)')'-rhess <filename>  read TM, ORCA or G09 hessian           '
  write(io,'(a)')'-hmass             mass-weight hessian for ANC'

src/getgrad.F90

@@ -185,12 +185,12 @@ subroutine getgrad
   aa='rm -f forcedump.dat force.dat > /dev/null'
   call system (trim(aa))
   call wamber(nat3,xyz,'amber.rst')
-  if(APBS) then
-  aaa='sander.APBS -O -i '//trim(ambin)// &
-               ' -c amber.rst'// &
-               ' -p '//trim(ambtop)// &
-               ' -o '//trim(xjob)
-  endif
+  ! if(APBS) then
+  ! aaa='sander.APBS -O -i '//trim(ambin)// &
+  !              ' -c amber.rst'// &
+  !              ' -p '//trim(ambtop)// &
+  !              ' -o '//trim(xjob)
+  ! endif
   call system(trim(command_amber))
   !call wamber(nat3,xyz)
   call ambergrad(nat,grad,energy,apbs)
@@ -813,15 +813,25 @@ subroutine drivergrad(xnat,grad,e)
 subroutine egradfile(xnat,grad,e,fname)
 use fiso, only: r8
 implicit none
-integer i,xnat,io
+integer i,xnat,io,c
 character(*) fname
 real(r8) grad(3,xnat),e
+logical da
 
+inquire(file=fname, exist=da)
+if(.not.da) then
+  call error('missing file: '//fname)
+  return
+endif
 open(newunit=io,file=fname)
-read(io,*) e
+read(io,*,end=666) e
+c=0
 do i=1,xnat
-   read(io,*) grad(1:3,i)
+  c=c+1
+   read(io,*,end=666) grad(1:3,i)
 enddo
+666 continue
+if (c/=xnat) call error('incomplete gradient')
 close(io)
 end subroutine
 

src/progs.f90

@@ -18,6 +18,8 @@ subroutine prepPROG
          write(stdout,*) " *  Gaussian                                                   *"
     elseif(amber) then
          write(stdout,*) " *  sander                                                     *"
+    elseif(gei) then
+         write(stdout,*) " *  custom interface                                  *"
     elseif(numgrad) then
            write(stdout,*) " * using parallel numerical gradients                             *"
          if(nproc>1) then