From 57b9d778900fd60a3a162dbc9ac9969f79c98a8c Mon Sep 17 00:00:00 2001 From: Vikrant Tripathy Date: Mon, 4 Nov 2024 01:11:41 -0800 Subject: [PATCH 1/7] enable the feature of skipping SCF by reading density matrix --- src/main.f90 | 3 ++- src/modules/quick_method_module.f90 | 12 +++++++++++- src/modules/quick_oeproperties_module.f90 | 21 +++++++++++++++++++- src/modules/quick_scf_module.f90 | 24 +++++++++++++++++------ src/subs/io.f90 | 18 +++++++++-------- 5 files changed, 61 insertions(+), 17 deletions(-) diff --git a/src/main.f90 b/src/main.f90 index c54faa945..00bdfe4a4 100644 --- a/src/main.f90 +++ b/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() diff --git a/src/modules/quick_method_module.f90 b/src/modules/quick_method_module.f90 index 7904ab0c0..87a5f8f3e 100644 --- a/src/modules/quick_method_module.f90 +++ b/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 TO 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")') @@ -655,7 +661,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 +902,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 diff --git a/src/modules/quick_oeproperties_module.f90 b/src/modules/quick_oeproperties_module.f90 index a438329d0..034b2badb 100644 --- a/src/modules/quick_oeproperties_module.f90 +++ b/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 diff --git a/src/modules/quick_scf_module.f90 b/src/modules/quick_scf_module.f90 index fd8761098..de73b345c 100644 --- a/src/modules/quick_scf_module.f90 +++ b/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) diff --git a/src/subs/io.f90 b/src/subs/io.f90 index 657f471de..e90a7a0eb 100644 --- a/src/subs/io.f90 +++ b/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 From 6279858d41a5ae8ff2c180cbf2aa971ad8bda1a2 Mon Sep 17 00:00:00 2001 From: Vikrant Tripathy Date: Mon, 4 Nov 2024 01:30:32 -0800 Subject: [PATCH 2/7] CUDA in ESP charge calculation --- src/modules/quick_oeproperties_module.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/modules/quick_oeproperties_module.f90 b/src/modules/quick_oeproperties_module.f90 index 034b2badb..02c1cb2cb 100644 --- a/src/modules/quick_oeproperties_module.f90 +++ b/src/modules/quick_oeproperties_module.f90 @@ -333,7 +333,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 @@ -343,7 +343,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. From 33394bcdda3251851515ef6192ed85f65ddbdaa3 Mon Sep 17 00:00:00 2001 From: Vikrant Tripathy Date: Fri, 8 Nov 2024 14:13:38 -0800 Subject: [PATCH 3/7] large no of grid points for ESP evaluation is now supported --- src/cuda/gpu_oeprop.h | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/cuda/gpu_oeprop.h b/src/cuda/gpu_oeprop.h index 3c445f09e..b4a03f701 100644 --- a/src/cuda/gpu_oeprop.h +++ b/src/cuda/gpu_oeprop.h @@ -192,13 +192,17 @@ __global__ void getOEPROP_kernel(){ unsigned int totalatom = devSim.natom+devSim.nextatom; unsigned int totalpoint = devSim.nextpoint; - for (QUICKULL i = offset; i < jshell * jshell * totalpoint; i+= totalThreads) { + QUICKULL jshellsq = (QUICKULL) (jshell * jshell); + long double inv_jshell2 = (long double) (1.0/jshellsq); + QUICKULL ncalcs = (QUICKULL) (jshellsq * totalpoint); + + for (QUICKULL i = offset; i < ncalcs; i+= totalThreads) { // use the global index to obtain shell pair. Note that here we obtain a couple of indices that helps us to obtain // shell number (ii and jj) and quantum numbers (iii, jjj). - - unsigned int ipoint = (int) i/(jshell * jshell); - unsigned int idx = i - ipoint * jshell * jshell; + + unsigned int ipoint = (unsigned int) (i*inv_jshell2); + unsigned int idx = (unsigned int) (i - ipoint * jshellsq); #ifdef CUDA_MPIV if(devSim.mpi_boeicompute[idx] > 0){ @@ -215,7 +219,6 @@ __global__ void getOEPROP_kernel(){ int iii = devSim.sorted_Qnumber[II]; int jjj = devSim.sorted_Qnumber[JJ]; - // compute coulomb attraction for the selected shell pair. iclass_oeprop(iii, jjj, ii, jj, ipoint, totalpoint, totalatom, devSim.YVerticalTemp+offset, devSim.store+offset, devSim.store2+offset); From b1d9b1fc8613630d85108892775361a46d247ff6 Mon Sep 17 00:00:00 2001 From: Vikrant Tripathy Date: Fri, 8 Nov 2024 14:13:58 -0800 Subject: [PATCH 4/7] long double replaced by double to remove warning --- src/cuda/gpu_oeprop.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cuda/gpu_oeprop.h b/src/cuda/gpu_oeprop.h index b4a03f701..257ccb9cf 100644 --- a/src/cuda/gpu_oeprop.h +++ b/src/cuda/gpu_oeprop.h @@ -193,7 +193,7 @@ __global__ void getOEPROP_kernel(){ unsigned int totalpoint = devSim.nextpoint; QUICKULL jshellsq = (QUICKULL) (jshell * jshell); - long double inv_jshell2 = (long double) (1.0/jshellsq); + double inv_jshell2 = (double) (1.0/jshellsq); QUICKULL ncalcs = (QUICKULL) (jshellsq * totalpoint); for (QUICKULL i = offset; i < ncalcs; i+= totalThreads) { From 91a780c8ce817bf3425f397208c7d27247119c4e Mon Sep 17 00:00:00 2001 From: Vikrant Tripathy Date: Fri, 8 Nov 2024 15:46:08 -0800 Subject: [PATCH 5/7] preposition error --- src/modules/quick_method_module.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modules/quick_method_module.f90 b/src/modules/quick_method_module.f90 index 87a5f8f3e..6320745aa 100644 --- a/src/modules/quick_method_module.f90 +++ b/src/modules/quick_method_module.f90 @@ -426,7 +426,7 @@ 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 TO 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")') From ed04cbfdc5ee9ba2c3e9cba9d7b4ff7cf327bf9f Mon Sep 17 00:00:00 2001 From: Vikrant Tripathy Date: Fri, 8 Nov 2024 17:52:59 -0800 Subject: [PATCH 6/7] QUICK can now read keywords from multiple lines --- src/read_job_and_atom.f90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/read_job_and_atom.f90 b/src/read_job_and_atom.f90 index ed3439d2f..c05f3f9ef 100644 --- a/src/read_job_and_atom.f90 +++ b/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 From b8e98762503b4830f9e6e6d4989311a01f8d0c00 Mon Sep 17 00:00:00 2001 From: Vikrant Tripathy Date: Sat, 9 Nov 2024 19:04:40 -0800 Subject: [PATCH 7/7] increased the max length of keyword string --- src/modules/quick_method_module.f90 | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/modules/quick_method_module.f90 b/src/modules/quick_method_module.f90 index 6320745aa..0159b002c 100644 --- a/src/modules/quick_method_module.f90 +++ b/src/modules/quick_method_module.f90 @@ -522,8 +522,6 @@ subroutine print_quick_method(self,io,ierr) end subroutine print_quick_method - - !------------------------ ! read quick_method !------------------------ @@ -532,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. @@ -1082,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 @@ -1123,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