Merge branch 'method_module' into oeproperties-dev

src/main.f90

@@ -268,8 +268,9 @@ program quick
                           -(timer_end%T2elb-timer_begin%T2elb)
 
     if (.not.quick_method%opt .and. .not.quick_method%grad) then
+      if(.not.quick_method%Skip)then
         SAFE_CALL(getEnergy(.false.,ierr))
-        
+      endif
         ! One electron properties (ESP, EField)
 
         !call generate_MKS_surfaces()

src/modules/quick_method_module.f90

@@ -44,7 +44,9 @@ module quick_method_module
         logical :: debug =  .false.    ! debug mode
         logical :: nodirect = .false.  ! conventional scf
         logical :: readDMX =  .false.  ! flag to read density matrix
+        logical :: readPMat = .false.  ! flag to read density matrix
         logical :: writePMat = .false. ! flag to write density matrix
+        logical :: Skip = .false.      ! Skip SCF and go to property calculations
         logical :: readSAD = .true.    ! flag to read SAD guess
         logical :: writeSAD = .false.  ! flag to write SAD guess
         logical :: diisSCF =  .false.  ! DIIS SCF
@@ -89,7 +91,7 @@ module quick_method_module
         integer :: iSG = 1             ! =0. SG0, =1. SG1(DEFAULT)
 
         ! Initial guess part
-        logical :: SAD = .true.        ! SAD initial guess(defualt
+        logical :: SAD = .true.        ! SAD initial guess(default)
         logical :: MFCC = .false.      ! MFCC
 
         ! this part is about ECP
@@ -246,7 +248,9 @@ subroutine broadcast_quick_method(self, ierr)
             call MPI_BCAST(self%calcDensLap,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%gridspacing,1,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%lapGridSpacing,1,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
+            call MPI_BCAST(self%readPMat,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%writePMat,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
+            call MPI_BCAST(self%Skip,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%extCharges,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%ext_grid,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%extgrid_angstrom,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
@@ -422,7 +426,9 @@ subroutine print_quick_method(self,io,ierr)
             if (self%printEnergy) write(io,'(" PRINT ENERGY EVERY CYCLE")')
 
             if (self%readDMX)   write(io,'(" READ DENSITY MATRIX FROM FILE")')
+            if (self%readPMat) write(io,'(" READ DENSITY MATRIX From FILE")')
             if (self%writePMat) write(io,'(" WRITE DENSITY MATRIX TO FILE")')
+            if (self%Skip) write(io,'(" SKIPPING SCF CALCULATION")')
             if (self%readSAD)   write(io,'(" READ SAD GUESS FROM FILE")')
             if (self%writeSAD)   write(io,'(" WRITE SAD GUESS TO FILE")')
 
@@ -516,8 +522,6 @@ subroutine print_quick_method(self,io,ierr)
 
         end subroutine print_quick_method
 
-
-
         !------------------------
         ! read quick_method
         !------------------------
@@ -526,13 +530,13 @@ subroutine read_quick_method(self,keywd,ierr)
             use quick_mpi_module
             use quick_files_module, only : write_molden
             implicit none
-            character(len=200) :: keyWD
-            character(len=200) :: tempstring
+            character(len=300) :: keyWD
+            character(len=300) :: tempstring
             integer :: itemp,i,j
             type (quick_method_type) self
             integer, intent(inout) :: ierr
 
-            call upcase(keyWD,200)
+            call upcase(keyWD,300)
             if (index(keyWD,'PDB').ne. 0)       self%PDB=.true.
             if (index(keyWD,'MFCC').ne.0)       self%MFCC=.true.
             if (index(keyWD,'FMM').ne.0)        self%FMM=.true.
@@ -655,7 +659,9 @@ subroutine read_quick_method(self,keywd,ierr)
             end if
             if (index(keyWD,'ZMAKE').ne.0)      self%zmat=.true.
             if (index(keyWD,'DIPOLE').ne.0)     self%dipole=.true.
+            if (index(keyWD,'READ').ne.0)      self%readPMat=.true.
             if (index(keyWD,'WRITE').ne.0)      self%writePMat=.true.
+            if (index(keyWD,'SKIP').ne.0)      self%Skip=.true.
             if (index(keyWD,'EXTCHARGES').ne.0) self%EXTCHARGES=.true.
             if (index(keyWD,'EXTGRID').ne.0) self%ext_grid=.true.
             !if (index(keyWD,'EXTGRID_ANGSTROM').ne.0) self%extgrid_angstrom=.true.
@@ -894,7 +900,9 @@ subroutine init_quick_method(self,ierr)
             self%hasF = .false.            ! If f orbitial is contained
             self%calcDens = .false.    ! calculate density
             self%calcDensLap = .false. ! calculate density lap
+            self%readPMat = .false.    ! Input density matrix
             self%writePMat = .false.   ! Output density matrix
+            self%Skip = .false.        ! Skipping single point SCF calculation
             self%extCharges = .false.  ! external charge
             self%ext_grid = .false.    ! external grid points
             self%extgrid_angstrom = .false.   ! external grid points (same as above) output in angstrom
@@ -1072,10 +1080,10 @@ subroutine set_libxc_func_info(f_keywd, self,ierr)
            use quick_exception_module
 
            implicit none
-           character(len=200), intent(in) :: f_keywd
+           character(len=300), intent(in) :: f_keywd
            type(quick_method_type), intent(inout) :: self
            integer, intent(inout) :: ierr
-           character(len=200) :: func1, func2 
+           character(len=300) :: func1, func2 
            character(len=256) :: functional_name
            integer :: f_id, nof_f, istart, iend, imid
            double precision :: x_hyb_coeff
@@ -1113,7 +1121,7 @@ subroutine set_libxc_func_info(f_keywd, self,ierr)
             .and. (index(functional_name,'mgga') .eq. 0))  then
                 functional_name=trim(functional_name)
 
-                call upcase(functional_name,200)
+                call upcase(functional_name,300)
 
                 if((trim(functional_name) == trim(func1)) .or. (trim(functional_name) == trim(func2))) then 
                         nof_f=nof_f+1

src/modules/quick_oeproperties_module.f90

@@ -29,15 +29,34 @@ module quick_oeproperties_module
 
  Subroutine compute_oeprop()
    use quick_method_module, only: quick_method
-   use quick_files_module, only : ioutfile
+   use quick_files_module, only : ioutfile, iDataFile, dataFileName
    use quick_molsurface_module, only: generate_MKS_surfaces
    use quick_molspec_module, only: quick_molspec
    use quick_mpi_module, only: master, mpierror
+   use quick_calculated_module, only: quick_qm_struct
 #ifdef MPIV
    use mpi
 #endif
    implicit none
 
+   logical fail
+   integer ierr, nbasis
+
+   if(quick_method%readPMat)then
+     nbasis = quick_molspec%nbasis
+     if(master)then
+       open(unit=iDataFile,file=dataFileName,status='OLD',form='UNFORMATTED')
+       call rchk_int(iDataFile, "nbasis", nbasis, fail)
+       call rchk_darray(iDataFile, "dense", nbasis, nbasis, 1, quick_qm_struct%dense, fail)
+       call rchk_darray(iDataFile, "denseSave", nbasis, nbasis, 1, quick_qm_struct%denseSave, fail)
+       close(iDataFile)
+     endif
+#ifdef MPIV
+     call MPI_BCAST(quick_qm_struct%denseSave,nbasis*nbasis,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
+     call MPI_BCAST(quick_qm_struct%dense,nbasis*nbasis,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
+#endif
+   endif
+
    if (quick_method%ext_grid) then
       call compute_oeprop_grid(quick_molspec%nextpoint,quick_molspec%extpointxyz)
    else if (quick_method%esp_charge) then
@@ -319,7 +338,7 @@ subroutine compute_ESP_charge(npoints,xyz_points,esp)
    allocate(IPIV(natom+1))
 
    LDA = natom+1
-   CALL dgetrf(natom+1,natom+1,A,LDA,IPIV,ierr)
+   CALL DGETRF(natom+1,natom+1,A,LDA,IPIV,ierr)
    call DGETRI(natom+1,A,LDA,IPIV,WORK,LWORK,ierr)
 
    if (ierr /= 0) then
@@ -329,7 +348,11 @@ subroutine compute_ESP_charge(npoints,xyz_points,esp)
 
 !  q = A-1*B
 
+#if defined CUDA
+   call CUBLAS_DGEMV('N',natom+1,natom+1,One,A,LDA,B,1,Zero,q,1)
+#else
    call DGEMV('N',natom+1,natom+1,One,A,LDA,B,1,Zero,q,1)
+#endif
 
 !  B is copied to charge array.
 

src/modules/quick_scf_module.f90

@@ -153,7 +153,7 @@ end subroutine scf
 
   ! electdiis
   !-------------------------------------------------------
-  ! 11/02/2010 Yipu Miao: Add paralle option for HF calculation
+  ! 11/02/2010 Yipu Miao: Add parallel option for HF calculation
   subroutine electdiis(jscf,ierr)
 
      use allmod
@@ -182,9 +182,11 @@ subroutine electdiis(jscf,ierr)
 #endif
 
      implicit none
-
+
+     logical fail
+
      ! variable inputed to return
-     integer :: jscf                ! scf interation
+     integer :: jscf                ! scf iteration
      integer, intent(inout) :: ierr
 
      logical :: diisdone = .false.  ! flag to indicate if diis is done
@@ -197,7 +199,7 @@ subroutine electdiis(jscf,ierr)
      double precision :: Sum2Mat,rms
      integer :: I,J,K,L,IERROR
 
-     double precision :: oldEnergy=0.0d0,E1e ! energy for last iteriation, and 1e-energy
+     double precision :: oldEnergy=0.0d0,E1e ! energy for last iteration, and 1e-energy
      double precision :: PRMS,PCHANGE, tmp
 
      double precision :: c_coords(3),c_zeta,c_chg
@@ -374,7 +376,7 @@ subroutine electdiis(jscf,ierr)
            !-----------------------------------------------
            ! The matrix multiplier comes from Steve Dixon. It calculates
            ! C = Transpose(A) B.  Thus to utilize this we have to make sure that the
-           ! A matrix is symetric. First, calculate DENSE*S and store in the scratch
+           ! A matrix is symmetric. First, calculate DENSE*S and store in the scratch
            ! matrix hold.Then calculate O*(DENSE*S).  As the operator matrix is symmetric, the
            ! above code can be used. Store this (the ODS term) in the all error
            ! matrix.
@@ -714,7 +716,17 @@ subroutine electdiis(jscf,ierr)
         !--------------- END MPI/ALL NODES -------------------------------------
 
         if (master) then
-
+
+          if(quick_method%writepmat)then
+             ! open data file then write calculated info to dat file
+             call quick_open(iDataFile, dataFileName, 'R', 'U', 'A',.true.,ierr)
+!             rewind(iDataFile)
+             call wchk_int(iDataFile, "nbasis", nbasis, fail)
+             call wchk_darray(iDataFile, "dense",    nbasis, nbasis, 1, quick_qm_struct%dense,    fail)
+             call wchk_darray(iDataFile, "denseSave",    nbasis, nbasis, 1, quick_qm_struct%denseSave,    fail)
+             close(iDataFile)
+          endif
+
 #ifdef USEDAT
            ! open data file then write calculated info to dat file
            SAFE_CALL(quick_open(iDataFile, dataFileName, 'R', 'U', 'R',.true.,ierr)

src/read_job_and_atom.f90

@@ -29,8 +29,8 @@ subroutine read_job_and_atom(ierr)
    ! Instant variables
    integer i,j,k
    integer istrt
-   character(len=200) :: keyWD
-   character(len=20) :: tempstring
+   character(len=300) :: keyWD
+   character(len=100) :: tempstring
    integer, intent(inout) :: ierr
 
    if (master) then
@@ -45,11 +45,16 @@ subroutine read_job_and_atom(ierr)
         keyWD=quick_api%Keywd
       else
         SAFE_CALL(quick_open(infile,inFileName,'O','F','W',.true.,ierr))
-        read (inFile,'(A200)') keyWD
+        keyWD(:)=''
+        do while(.true.)
+          read (inFile,'(A100)') tempstring
+          if(trim(tempstring).eq.'') exit
+          keyWD=trim(keyWD)//' '//trim(tempstring)
+        end do
       endif
 
-      call upcase(keyWD,200)
-      write(iOutFile,'("  KEYWORD=",a)') keyWD
+      call upcase(keyWD,300)
+      write(iOutFile,'("  KEYWORD=",a)') trim(keyWD)
 
       ! These interfaces,"read","check" and "print" are from quick_method_module
       SAFE_CALL(read(quick_method,keyWD,ierr))     ! read method from Keyword

src/subs/io.f90

@@ -35,6 +35,7 @@ subroutine wchk_int(chk,key,nvalu,fail)
    write(chk) 'I '
    write(chk) nvalu
    fail=1
+   return
    200  return
 
 end
@@ -376,18 +377,19 @@ subroutine wchk_darray(chk,key,x,y,z,dim,fail)
       enddo
    endif
 
-   fail=0
-   do
-      read(chk,end=100,err=200)
-   enddo
-
-   100  rewind(chk)
+!   fail=0
+!   do
+!      read(chk,end=100,err=200)
+!   enddo
+!
+!   100  rewind(chk)
    write(chk) '#'//kline(1:40)
    write(chk) 'RR'
    write(chk) x*y*z
    write(chk) (((dim(i,j,k),i=1,x),j=1,y),k=1,z)
-   fail=1
-   200  return
+!   fail=1
+!   200  return
+   return
 
 end