Initital working implementation / serial version

src/CMakeLists.txt

@@ -33,7 +33,7 @@ set(QUICK_MODULES_SOURCES
 	quick_uscf_module.f90 quick_sad_guess_module.f90 quick_optimizer_module.f90
         quick_cew_module.f90 quick_lri_module.f90 quick_lri_grad_module.f90
 	oshell_quick_eri_module.f90 oshell_quick_eri_grad_module.f90
-	oshell_quick_gradient_module.f90 quick_dftd3_module.f90 quick_molden_module.f90)
+	oshell_quick_gradient_module.f90 quick_dftd3_module.f90 quick_molden_module.f90 quick_oeproperties_module.f90)
 
 set(QUICK_SUBS_SOURCES Angles.f90 copyDMat.f90 copySym.f90
 	degen.f90 denspt.f90 diag.f90 dipole.f90

src/main.f90

@@ -30,6 +30,7 @@ program quick
     use quick_exception_module
     use quick_cshell_eri_module, only: getEriPrecomputables
     use quick_cshell_gradient_module, only: cshell_gradient
+    use quick_oeproperties_module, only: compute_esp
     use quick_oshell_gradient_module, only: oshell_gradient
     use quick_optimizer_module
     use quick_sad_guess_module, only: getSadGuess
@@ -39,6 +40,8 @@ program quick
     use mpi
 #endif
 
+
+
     implicit none
 
 #if defined CUDA || defined HIP
@@ -52,6 +55,9 @@ program quick
     integer :: i,j,k
     double precision t1_t, t2_t
     common /timer/ t1_t, t2_t
+    double precision, dimension(:), allocatable :: esp_array
+
+
     !------------------------------------------------------------------
     ! 1. The first thing that must be done is to initialize and prepare files
     !------------------------------------------------------------------
@@ -225,6 +231,8 @@ program quick
 
       call gpu_upload_oei(quick_molspec%nExtAtom, quick_molspec%extxyz, quick_molspec%extchg, ierr)
 
+      call gpu_upload_oei(quick_molspec%nextpoint, quick_molspec%extxyz, ierr)
+
     endif
 #endif
 
@@ -286,8 +294,6 @@ program quick
        end if
     endif
 
-
-
     !------------------------------------------------------------------
     ! 6. Other job options
     !-----------------------------------------------------------------
@@ -344,6 +350,11 @@ program quick
 
     endif
 
+    ! 6.e Electrostatic Potential
+    if (quick_method%esp_grid) then
+        call compute_esp(ierr)
+    end if
+
     ! Now at this point we have an energy and a geometry.  If this is
     ! an optimization job, we now have the optimized geometry.
 

src/modules/Makefile

@@ -37,7 +37,7 @@ modobj=$(objfolder)/quick_exception_module.o $(objfolder)/quick_mpi_module.o $(o
 		$(objfolder)/quick_eri_module.o $(objfolder)/oshell_quick_eri_module.o $(objfolder)/quick_lri_module.o \
 		$(objfolder)/quick_molden_module.o $(objfolder)/quick_eri_grad_module.o \
 		$(objfolder)/oshell_quick_eri_grad_module.o $(objfolder)/quick_lri_grad_module.o \
-		$(objfolder)/quick_cew_module.o $(objfolder)/quick_oei_module.o $(objfolder)/quick_dft_module.o \
+		$(objfolder)/quick_cew_module.o $(objfolder)/quick_oei_module.o $(objfolder)/quick_oeproperties_module.o $(objfolder)/quick_dft_module.o \
 		$(objfolder)/quick_scf_operator_module.o $(objfolder)/quick_scf_module.o $(objfolder)/quick_gradient_module.o \
 		$(objfolder)/oshell_quick_gradient_module.o \
 		$(objfolder)/quick_uscf_operator_module.o $(objfolder)/quick_uscf_module.o $(objfolder)/quick_sad_guess_module.o \

src/modules/quick_calculated_module.f90

@@ -167,6 +167,9 @@ module quick_calculated_module
       ! Mulliken charge and Lowdin charge
       double precision,dimension(:),allocatable :: Mulliken,Lowdin
 
+      ! Electrostatic potential 
+      double precision,dimension(:), allocatable   :: ESP
+
 
    end type quick_qm_struct_type
 
@@ -209,9 +212,9 @@ module quick_calculated_module
    !----------------------
 contains
 
-   !--------------
+   !--------------------------------------
    ! subroutine to allocate variables
-   !--------------
+   !---------------------------------------
    subroutine allocate_quick_qm_struct(self)
       use quick_method_module,only: quick_method
       use quick_molspec_module,only: quick_molspec

src/modules/quick_files_module.f90

@@ -37,6 +37,8 @@ module quick_files_module
     character(len=80) :: dataFileName   = ''
     character(len=80) :: intFileName    = ''
     character(len=80) :: moldenFileName = ''
+    character(len=80) :: propFileName   = ''
+
 
     ! Basis set and directory
     character(len=240) :: basisDir       = ''
@@ -64,13 +66,18 @@ module quick_files_module
     integer :: iDataFile      = DATAFILEHANDLE   ! Data file, similar to chk file in gaussian
     integer :: iIntFile       = INTFILEHANDLE    ! integral file
     integer :: iMoldenFile    = MOLDENFILEHANDLE ! molden file
+    integer :: iPropFile    = PROPFILEHANDLE     ! properties file
 
     logical :: fexist = .false.         ! Check if file exists
 
     logical :: isTemplate = .false.   ! is input file a template (i.e. only the keywords)
     integer :: wrtStep = 1            ! current step for writing to output file.
     logical :: write_molden = .false. ! flag to export data into molden format
 
+    ! is this necessary?
+   ! logical :: write_prop = .false. ! flag to export data into prop format
+
+
     contains
 
     !------------
@@ -123,6 +130,7 @@ subroutine set_quick_files(api,ierr)
         dataFileName=inFileName(1:i-1)//'.dat'
         intFileName=inFileName(1:i-1)//'.int'
         moldenFileName=inFileName(1:i-1)//'.molden'
+        propFileName=inFileName(1:i-1)//'.prop'
 
         return
 
@@ -279,5 +287,4 @@ subroutine print_quick_io_file(io,ierr)
     end subroutine print_quick_io_file
 
 
-
 end module quick_files_module

src/modules/quick_method_module.f90

@@ -53,6 +53,9 @@ module quick_method_module
         logical :: grad = .false.      ! if calculate gradient
         logical :: analGrad =  .false. ! Analytical Gradient
         logical :: analHess =  .false. ! Analytical Hessian Matrix
+        logical :: esp_grid =  .false.      ! Electrostatic potential on a grid (ESP)
+        logical :: efield_grid =  .false.   ! Electrostatic field (EFIELD)
+        logical :: efg_grid =  .false.      ! Electrostatic field gradient (EFG)
         logical :: diisOpt =  .false.  ! DIIS Optimization
         logical :: core =  .false.     ! Add core
         logical :: annil =  .false.    ! Annil Spin Contamination
@@ -70,7 +73,8 @@ module quick_method_module
                                        ! Density lapcacian file gridspacing
 
         logical :: PDB = .false.       ! PDB input
-        logical :: extCharges = .false.! external charge
+        logical :: extCharges = .false.! external charge (x,y,z,q)
+        logical :: ext_grid = .false.  ! external grid points (x,y,z)
 
 
         ! those methods are mostly for research use
@@ -98,7 +102,7 @@ module quick_method_module
         integer :: maxdiisscf = 10
 
         ! start cycle for delta density cycle
-        integer :: ncyc =3
+        integer :: ncyc = 3
 
         ! following are some cutoff criteria
         double precision :: coreIntegralCutoff = 1.0d-12 ! cutoff for 1e integral prescreening
@@ -219,6 +223,9 @@ subroutine broadcast_quick_method(self, ierr)
             call MPI_BCAST(self%grad,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%analGrad,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%analHess,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
+            call MPI_BCAST(self%esp_grid,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
+            call MPI_BCAST(self%efield_grid,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
+            call MPI_BCAST(self%efg_grid,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%diisOpt,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%core,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%annil,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
@@ -236,6 +243,7 @@ subroutine broadcast_quick_method(self, ierr)
             call MPI_BCAST(self%lapGridSpacing,1,mpi_double_precision,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%writePMat,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%PDB,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%SAD,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
             call MPI_BCAST(self%FMM,1,mpi_logical,0,MPI_COMM_WORLD,mpierror)
@@ -300,7 +308,7 @@ subroutine print_quick_method(self,io,ierr)
             if (io.ne.0) then
             write(io,'(" ============== JOB CARD =============")')
             if (self%HF) then
-                write(io,'(" METHOD = HATREE FOCK")')
+                write(io,'(" METHOD = HARTREE FOCK")')
             else if (self%MP2) then
                 write(io,'(" METHOD = SECOND ORDER PERTURBATION THEORY")')
             else if (self%DFT) then
@@ -418,6 +426,7 @@ subroutine print_quick_method(self,io,ierr)
             if (self%custECP)   write(io,'(" CUSTOM ECP BASIS SET")')
 
             if (self%extCharges)write(io,'(" EXTERNAL CHARGES")')
+            if (self%ext_grid)write(io,'(" EXTERNAL GRID POINTS")')
             if (self%core)      write(io,'(" SUM INNER ELECTRONS INTO CORE")')
             if (self%debug)     write(io,'(" DEBUG MODE")')
 
@@ -454,6 +463,11 @@ subroutine print_quick_method(self,io,ierr)
 
             if (self%grad)      write(io,'(" GRADIENT CALCULATION")')
 
+            ! computing esp, efield and efg
+            if (self%esp_grid)      write(io,'(" ELECTROSTATIC POTENTIAL CALCULATION")')
+            if (self%efield_grid)      write(io,'(" ELECTROSTATIC FIELD CALCULATION")')
+            if (self%efg_grid)      write(io,'(" ELECTROSTATIC FIELD GRADIENT CALCULATION")')
+
             if (self%DIVCON) then
                 write(io,'(" DIV & CON METHOD")',advance="no")
                 if (self%ifragbasis .eq. 1) then
@@ -632,6 +646,7 @@ subroutine read_quick_method(self,keywd,ierr)
             if (index(keyWD,'DIPOLE').ne.0)     self%dipole=.true.
             if (index(keyWD,'WRITE').ne.0)      self%writePMat=.true.
             if (index(keyWD,'EXTCHARGES').ne.0) self%EXTCHARGES=.true.
+            if (index(keyWD,'EXTGRID').ne.0) self%ext_grid=.true.
             if (index(keyWD,'FORCE').ne.0)      self%grad=.true.
 
             if (index(keyWD,'NODIRECT').ne.0)      self%NODIRECT=.true.
@@ -780,6 +795,25 @@ subroutine read_quick_method(self,keywd,ierr)
                     ierr=35
                 endif
             endif
+
+
+          ! Electrostatic Potential(ESP), Field (EF), Field Gradient (EFG) on a grid
+           if (index(keyWD,'EFG').ne.0) then
+               self%esp_grid=.true.
+               self%efield_grid=.true.
+               self%efg_grid=.true.
+               self%ext_grid=.true.
+           endif
+           if (index(keyWD,'EFIELD').ne.0) then
+               self%esp_grid=.true.
+               self%efield_grid=.true.
+               self%ext_grid=.true.
+           endif
+           if (index(keyWD,'ESP').ne.0) then
+               self%esp_grid=.true.
+               self%ext_grid=.true.
+           endif
+
             CHECK_ERROR(ierr)
 
         end subroutine read_quick_method
@@ -813,17 +847,20 @@ subroutine init_quick_method(self,ierr)
             self%PBSOL = .false.     ! PB Solvent
             self%UNRST =  .false.    ! Unrestricted
 
-            self%debug =  .false.    ! debug mode
-            self%nodirect = .false.  ! conventional SCF
-            self%readDMX =  .false.  ! flag to read density matrix
-            self%readSAD =  .true.   ! flag to read sad guess
-            self%writeSAD = .false.  ! flag to write sad guess
-            self%diisSCF =  .false.  ! DIIS SCF
-            self%prtGap =  .false.   ! flag to print HOMO-LUMO gap
-            self%opt =  .false.      ! optimization
-            self%grad =  .false.     ! gradient
-            self%analGrad =  .true. ! Analytical Gradient
-            self%analHess =  .false. ! Analytical Hessian Matrix
+            self%debug =  .false.     ! debug mode
+            self%nodirect = .false.   ! conventional SCF
+            self%readDMX =  .false.   ! flag to read density matrix
+            self%readSAD =  .true.    ! flag to read sad guess
+            self%writeSAD = .false.   ! flag to write sad guess
+            self%diisSCF =  .false.   ! DIIS SCF
+            self%prtGap =  .false.    ! flag to print HOMO-LUMO gap
+            self%opt =  .false.       ! optimization
+            self%grad =  .false.      ! gradient
+            self%analGrad =  .true.   ! Analytical Gradient
+            self%analHess =  .false.  ! Analytical Hessian Matrix
+            self%esp_grid = .false.        ! Electrostatic potential (ESP) on grid
+            self%efield_grid = .false.     ! Electric field (EF) corresponding to ESP 
+            self%efg_grid = .false.        ! Electric field gradient (EFG)
 
             self%diisOpt =  .false.  ! DIIS Optimization
             self%core =  .false.     !
@@ -839,6 +876,7 @@ subroutine init_quick_method(self,ierr)
             self%calcDensLap = .false. ! calculate density lap
             self%writePMat = .false.   ! Output density matrix
             self%extCharges = .false.  ! external charge
+            self%ext_grid = .false.    ! external grid points
             self%PDB = .false.         ! PDB input
             self%SAD = .true.          ! SAD initial guess
             self%FMM = .false.         ! Fast Multipole