From c86f256192fc4860f043cd0fdec72cb8959b136f Mon Sep 17 00:00:00 2001 From: "Igor S. Gerasimov" <> Date: Fri, 13 May 2022 16:47:26 +0900 Subject: [PATCH] Remove defensive compilation --- Makefile | 6 +- O1.f90 | 169 ---------------------------------------------- otherfunc2.f90 | 180 ++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 179 insertions(+), 176 deletions(-) delete mode 100644 O1.f90 diff --git a/Makefile b/Makefile index a940282..c11adb7 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,5 @@ SIMD = -msse3 OPT = -O2 -qopenmp -qopenmp-link=static -threads -qopt-matmul $(SIMD) -diag-disable 8290,8291,6371,10316 -fpp -mkl -static-intel -OPT1 = -O1 -qopenmp -qopenmp-link=static -threads $(SIMD) -diag-disable 8290,8291,6371,10316 -fpp -mkl -static-intel #Options in the next line is for debugging purpose #OPTDBG = -O0 -qopenmp -diag-disable 8290,8291,6371 -threads -qopenmp-link=static -debug all -g -traceback -check all -fstack-protector -fpp -mkl -static-intel @@ -15,7 +14,7 @@ LIBRETAPATH = ./libreta_hybrid objects = define.o util.o plot.o Bspline.o sym.o libreta.o function.o GUI.o sub.o integral.o Lebedev-Laikov.o \ DFTxclib.o edflib.o fparser.o fileIO.o spectrum.o DOS.o Multiwfn.o 0123dim.o LSB.o \ population.o orbcomp.o bondorder.o topology.o excittrans.o otherfunc.o \ -otherfunc2.o otherfunc3.o O1.o surfana.o procgriddata.o AdNDP.o fuzzy.o CDA.o basin.o \ +otherfunc2.o otherfunc3.o surfana.o procgriddata.o AdNDP.o fuzzy.o CDA.o basin.o \ orbloc.o visweak.o EDA.o CDFT.o ETS_NOCV.o atmraddens.o NAONBO.o grid.o PBC.o hyper_polar.o deloc_aromat.o \ minpack.o blockhrr_012345.o ean.o hrr_012345.o eanvrr_012345.o boysfunc.o naiveeri.o ryspoly.o \ inquire.o @@ -153,9 +152,6 @@ otherfunc2.o : otherfunc2.f90 $(modules) otherfunc3.o : otherfunc3.f90 $(modules) $(FC) $(OPT) -c otherfunc3.f90 - -O1.o : O1.f90 $(modules) - $(FC) $(OPT1) -c O1.f90 surfana.o : surfana.f90 $(modules) $(FC) $(OPT) -c surfana.f90 diff --git a/O1.f90 b/O1.f90 deleted file mode 100644 index 1dee462..0000000 --- a/O1.f90 +++ /dev/null @@ -1,169 +0,0 @@ -!################################################################################## -!Some codes do not work normally if compiling via ifort -O2 in Linux possibly due to -!bug in compiler. These codes are collected in this file, they will be compiled via -O1 in Linux. -!In Windows, this file is also compiled with /O2 like other files. -!################################################################################## - -!!----------- Calculate center, first/second moments, radius of gyration, and
of a function -!Via ifort 2018 and 2021 with -O2, result obtained by option 1 is zero or NaN -subroutine funcmoment - use defvar - use util - use function - implicit none - real*8 intval, moment1(3), moment2(3, 3), moment2nuc(3, 3), funcval(radpot*sphpot), beckeweigrid(radpot*sphpot) - real*8 eigvecmat(3, 3), eigval(3) - type(content) gridatm(radpot*sphpot), gridatmorg(radpot*sphpot) - character selectyn - real*8 cenx, ceny, cenz, realcenx, realceny, realcenz - real*8 tmpval, xtmp, ytmp, ztmp - integer ifunc, isel, iradcut - integer i, iatm - integer iwalltime1, iwalltime2 - ifunc = 1 - cenx = 0 - ceny = 0 - cenz = 0 - do while (.true.) - write (*, *) - write (*, *) "0 Return" - write (*, *) "1 Calculate various quantities of the selected function" - write (*, *) "2 Calculate center and integral of the selected function" - write (*, "(a,i5)") " 3 Select the function to be studied, current:", ifunc - write (*, "(a,3f11.5,' Ang')") " 4 Set the center for option 1, current:", cenx*b2a, ceny*b2a, cenz*b2a - write (*, *) "5 Calculate center and integral of the absolute of the selected function" - read (*, *) isel - - if (isel == 0) then - return - else if (isel == 3) then - call selfunc_interface(1, ifunc) - cycle - else if (isel == 4) then - write (*, *) "Input X,Y,Z of the center in Angstrom, e.g. 2.0,0,1.5" - read (*, *) cenx, ceny, cenz - cenx = cenx/b2a !To Bohr - ceny = ceny/b2a - cenz = cenz/b2a - cycle - end if - - write (*, "(' Radial points:',i5,' Angular points:',i5,' Total:',i10,' per center')") radpot, sphpot, radpot*sphpot - call gen1cintgrid(gridatmorg, iradcut) - - call walltime(iwalltime1) - - intval = 0 - moment1 = 0 - moment2 = 0 - realcenx = 0 - realceny = 0 - realcenz = 0 - do iatm = 1, ncenter - write (*, "(' Processing center',i6,'(',a2,') /',i6)") iatm, a(iatm)%name, ncenter - gridatm%x = gridatmorg%x + a(iatm)%x !Move quadrature point to actual position in molecule - gridatm%y = gridatmorg%y + a(iatm)%y - gridatm%z = gridatmorg%z + a(iatm)%z -!$OMP parallel do shared(funcval) private(i) num_threads(nthreads) - do i = 1 + iradcut*sphpot, radpot*sphpot - funcval(i) = calcfuncall(ifunc, gridatm(i)%x, gridatm(i)%y, gridatm(i)%z) - end do -!$OMP end parallel do - call gen1cbeckewei(iatm, iradcut, gridatm, beckeweigrid, covr_tianlu, 3) - do i = 1 + iradcut*sphpot, radpot*sphpot - tmpval = funcval(i)*gridatmorg(i)%value*beckeweigrid(i) - xtmp = gridatm(i)%x - cenx - ytmp = gridatm(i)%y - ceny - ztmp = gridatm(i)%z - cenz - if (isel == 5) tmpval = abs(tmpval) - intval = intval + tmpval - moment1(1) = moment1(1) + xtmp*tmpval - moment1(2) = moment1(2) + ytmp*tmpval - moment1(3) = moment1(3) + ztmp*tmpval - moment2(1, 1) = moment2(1, 1) + xtmp*xtmp*tmpval - moment2(2, 2) = moment2(2, 2) + ytmp*ytmp*tmpval - moment2(3, 3) = moment2(3, 3) + ztmp*ztmp*tmpval - moment2(1, 2) = moment2(1, 2) + xtmp*ytmp*tmpval - moment2(2, 3) = moment2(2, 3) + ytmp*ztmp*tmpval - moment2(1, 3) = moment2(1, 3) + xtmp*ztmp*tmpval - realcenx = realcenx + gridatm(i)%x*tmpval - realceny = realceny + gridatm(i)%y*tmpval - realcenz = realcenz + gridatm(i)%z*tmpval - end do - end do - call walltime(iwalltime2) - write (*, "(' Calculation took up wall clock time',i10,' s',/)") iwalltime2 - iwalltime1 - - if (isel == 1) then - moment2(3, 1) = moment2(1, 3) - moment2(2, 1) = moment2(1, 2) - moment2(3, 2) = moment2(2, 3) - write (*, *) "Note: All data shown below are in a.u." - write (*, "(/,' Integral over whole space:',1PE16.8,/)") intval - write (*, "(' The first moment:')") - write (*, "(' X= ',1PE16.8,' Y= ',1PE16.8,' Z= ',1PE16.8)") moment1 - write (*, "(' Norm= ',1PE16.8,/)") sum(moment1**2) - write (*, "(' The second moment:')") - write (*, "(' XX=',1PE16.8,' XY=',1PE16.8,' XZ=',1PE16.8)") moment2(1, :) - write (*, "(' YX=',1PE16.8,' YY=',1PE16.8,' YZ=',1PE16.8)") moment2(2, :) - write (*, "(' ZX=',1PE16.8,' ZY=',1PE16.8,' ZZ=',1PE16.8)") moment2(3, :) - - call diagmat(moment2, eigvecmat, eigval, 300, 1D-10) - call sort(eigval) - write (*, "(a,3(1PE16.8))") ' Eigenvalues:', eigval - write (*, "(a,1PE16.8)") " Sum of eigenvalues (trace of the second moment tensor):", sum(eigval) - write (*, "(' Anisotropy:',1PE16.8,/)") eigval(3) - (eigval(1) + eigval(2))/2D0 - write (*, "(' Radius of gyration:',1PE16.8)") dsqrt((moment2(1, 1) + moment2(2, 2) + moment2(3, 3))/intval) - write (*, "(/,a,f16.6)") " Spatial extent of the function
:", sum(eigval) - - if (ifunc == 1) then - moment2nuc = 0 - do iatm = 1, ncenter - xtmp = a(iatm)%x - cenx - ytmp = a(iatm)%y - ceny - ztmp = a(iatm)%z - cenz - tmpval = a(iatm)%charge - moment2nuc(1, 1) = moment2nuc(1, 1) + xtmp*xtmp*tmpval - moment2nuc(2, 2) = moment2nuc(2, 2) + ytmp*ytmp*tmpval - moment2nuc(3, 3) = moment2nuc(3, 3) + ztmp*ztmp*tmpval - moment2nuc(1, 2) = moment2nuc(1, 2) + xtmp*ytmp*tmpval - moment2nuc(2, 3) = moment2nuc(2, 3) + ytmp*ztmp*tmpval - moment2nuc(1, 3) = moment2nuc(1, 3) + xtmp*ztmp*tmpval - end do - moment2nuc(3, 1) = moment2nuc(1, 3) - moment2nuc(2, 1) = moment2nuc(1, 2) - moment2nuc(3, 2) = moment2nuc(2, 3) - write (*, *) - write (*, "(' The quadrupole moment of nuclear charges:')") - write (*, "(' XX=',f16.8,' XY=',f16.8,' XZ=',f16.8)") moment2nuc(1, :) - write (*, "(' YX=',f16.8,' YY=',f16.8,' YZ=',f16.8)") moment2nuc(2, :) - write (*, "(' ZX=',f16.8,' ZY=',f16.8,' ZZ=',f16.8)") moment2nuc(3, :) - write (*, *) - write (*, "(' The quadrupole moment of the system:')") - write (*, "(' XX=',f16.8,' XY=',f16.8,' XZ=',f16.8)") moment2nuc(1, :) - moment2(1, :) - write (*, "(' YX=',f16.8,' YY=',f16.8,' YZ=',f16.8)") moment2nuc(2, :) - moment2(2, :) - write (*, "(' ZX=',f16.8,' ZY=',f16.8,' ZZ=',f16.8)") moment2nuc(3, :) - moment2(3, :) - end if - - else if (isel == 2 .or. isel == 5) then - realcenx = realcenx/intval - realceny = realceny/intval - realcenz = realcenz/intval - if (isel == 2) then - write (*, "(' Integral of the function:',1PE16.8,' a.u.')") intval - write (*, "(/,' Center of the function:')") - else - write (*, "(' Integral of the absolute of the function:',1PE16.8,' a.u.')") intval - write (*, "(/,' Center of the absolute of the function:')") - end if - write (*, "(' X=',f16.8,' Y=',f16.8,' Z=',f16.8,' Angstrom',/)") realcenx*b2a, realceny*b2a, realcenz*b2a - write (*, *) "Use this center for subsequent calculations? (y/n)" - read (*, *) selectyn - if (selectyn == 'y') then - cenx = realcenx - ceny = realceny - cenz = realcenz - end if - end if - end do -end subroutine diff --git a/otherfunc2.f90 b/otherfunc2.f90 index a99adb1..8e5fa79 100644 --- a/otherfunc2.f90 +++ b/otherfunc2.f90 @@ -1886,7 +1886,7 @@ end subroutine -!!--------- Calculate core-valence bifurcation (CVB) index and related quantities +!!--------- Calculate core-valence bifurcation (CVB) index and related quantities subroutine CVB_index use defvar use function @@ -1896,7 +1896,7 @@ real*8 ELF_x(nptELFcurve),ELF_y(nptELFcurve) write(*,*) "Original paper of CVB index: Theor. Chem. Acc., 104, 13 (2000)" write(*,*) -write(*,*) "------ Calculating core-valence bifurcation (CVB) and related quantities -----" +write(*,*) "------ Calculating core-valence bifurcation (CVB) and related quantities -----" write(*,*) "Input index of donor atom, hydrogen and acceptor atom in the H-bond (D-H...A)" write(*,*) "For example: 1,3,4" read(*,*) iD,iH,iA @@ -2895,3 +2895,179 @@ else end if end subroutine + +!!----------- Calculate center, first/second moments, radius of gyration, and
of a function +subroutine funcmoment + use defvar, only: nthreads + use defvar, only: content + use defvar, only: a, ncenter, b2a + use defvar, only: radpot, sphpot + use defvar, only: covr_tianlu + use util, only: walltime, diagmat, sort + use function, only: calcfuncall + implicit none + real*8 intval, moment1(3), moment2(3, 3), moment2nuc(3, 3) + real*8, allocatable :: funcval(:), beckeweigrid(:) + real*8 eigvecmat(3, 3), eigval(3) + type(content), allocatable :: gridatm(:), gridatmorg(:) + character selectyn + real*8 cenx, ceny, cenz, realcenx, realceny, realcenz + real*8 tmpval, xtmp, ytmp, ztmp + integer ifunc, isel, iradcut + integer i, iatm + integer iwalltime1, iwalltime2 + allocate(gridatm(radpot*sphpot), gridatmorg(radpot*sphpot)) + allocate(funcval(radpot*sphpot), beckeweigrid(radpot*sphpot)) + ifunc = 1 + cenx = 0 + ceny = 0 + cenz = 0 + do while (.true.) + write (*, *) + write (*, *) "0 Return" + write (*, *) "1 Calculate various quantities of the selected function" + write (*, *) "2 Calculate center and integral of the selected function" + write (*, "(a,i5)") " 3 Select the function to be studied, current:", ifunc + write (*, "(a,3f11.5,' Ang')") " 4 Set the center for option 1, current:", cenx*b2a, ceny*b2a, cenz*b2a + write (*, *) "5 Calculate center and integral of the absolute of the selected function" + read (*, *) isel + + if (isel == 0) then + return + else if (isel == 3) then + call selfunc_interface(1, ifunc) + cycle + else if (isel == 4) then + write (*, *) "Input X,Y,Z of the center in Angstrom, e.g. 2.0,0,1.5" + read (*, *) cenx, ceny, cenz + cenx = cenx/b2a !To Bohr + ceny = ceny/b2a + cenz = cenz/b2a + cycle + end if + + write (*, "(' Radial points:',i5,' Angular points:',i5,' Total:',i10,' per center')") radpot, sphpot, radpot*sphpot + call gen1cintgrid(gridatmorg, iradcut) + + call walltime(iwalltime1) + + intval = 0 + moment1 = 0 + moment2 = 0 + realcenx = 0 + realceny = 0 + realcenz = 0 + do iatm = 1, ncenter + write (*, "(' Processing center',i6,'(',a2,') /',i6)") iatm, a(iatm)%name, ncenter + gridatm%x = gridatmorg%x + a(iatm)%x !Move quadrature point to actual position in molecule + gridatm%y = gridatmorg%y + a(iatm)%y + gridatm%z = gridatmorg%z + a(iatm)%z +!$OMP parallel do shared(funcval) private(i) num_threads(nthreads) + do i = 1 + iradcut*sphpot, radpot*sphpot + funcval(i) = calcfuncall(ifunc, gridatm(i)%x, gridatm(i)%y, gridatm(i)%z) + end do +!$OMP end parallel do + call gen1cbeckewei(iatm, iradcut, gridatm, beckeweigrid, covr_tianlu, 3) +!$OMP parallel do num_threads(nthreads) & +!$OMP default(none) & +!$OMP private(i, xtmp, ytmp, ztmp, tmpval) & +!$OMP shared(iradcut, radpot, sphpot, isel, funcval, gridatm, gridatmorg, beckeweigrid, cenx, ceny, cenz) & +!$OMP reduction(+:intval, moment1, moment2, realcenx, realceny, realcenz) + do i = 1 + iradcut*sphpot, radpot*sphpot + tmpval = funcval(i)*gridatmorg(i)%value*beckeweigrid(i) + xtmp = gridatm(i)%x - cenx + ytmp = gridatm(i)%y - ceny + ztmp = gridatm(i)%z - cenz + if (isel == 5) tmpval = abs(tmpval) + intval = intval + tmpval + moment1(1) = moment1(1) + xtmp*tmpval + moment1(2) = moment1(2) + ytmp*tmpval + moment1(3) = moment1(3) + ztmp*tmpval + moment2(1, 1) = moment2(1, 1) + xtmp*xtmp*tmpval + moment2(2, 2) = moment2(2, 2) + ytmp*ytmp*tmpval + moment2(3, 3) = moment2(3, 3) + ztmp*ztmp*tmpval + moment2(1, 2) = moment2(1, 2) + xtmp*ytmp*tmpval + moment2(2, 3) = moment2(2, 3) + ytmp*ztmp*tmpval + moment2(1, 3) = moment2(1, 3) + xtmp*ztmp*tmpval + realcenx = realcenx + gridatm(i)%x*tmpval + realceny = realceny + gridatm(i)%y*tmpval + realcenz = realcenz + gridatm(i)%z*tmpval + end do +!$OMP end parallel do + end do + call walltime(iwalltime2) + write (*, "(' Calculation took up wall clock time',i10,' s',/)") iwalltime2 - iwalltime1 + + if (isel == 1) then + moment2(3, 1) = moment2(1, 3) + moment2(2, 1) = moment2(1, 2) + moment2(3, 2) = moment2(2, 3) + write (*, *) "Note: All data shown below are in a.u." + write (*, "(/,' Integral over whole space:',1PE16.8,/)") intval + write (*, "(' The first moment:')") + write (*, "(' X= ',1PE16.8,' Y= ',1PE16.8,' Z= ',1PE16.8)") moment1 + write (*, "(' Norm= ',1PE16.8,/)") sum(moment1**2) + write (*, "(' The second moment:')") + write (*, "(' XX=',1PE16.8,' XY=',1PE16.8,' XZ=',1PE16.8)") moment2(1, :) + write (*, "(' YX=',1PE16.8,' YY=',1PE16.8,' YZ=',1PE16.8)") moment2(2, :) + write (*, "(' ZX=',1PE16.8,' ZY=',1PE16.8,' ZZ=',1PE16.8)") moment2(3, :) + + call diagmat(moment2, eigvecmat, eigval, 300, 1D-10) + call sort(eigval) + write (*, "(a,3(1PE16.8))") ' Eigenvalues:', eigval + write (*, "(a,1PE16.8)") " Sum of eigenvalues (trace of the second moment tensor):", sum(eigval) + write (*, "(' Anisotropy:',1PE16.8,/)") eigval(3) - (eigval(1) + eigval(2))/2D0 + write (*, "(' Radius of gyration:',1PE16.8)") dsqrt((moment2(1, 1) + moment2(2, 2) + moment2(3, 3))/intval) + write (*, "(/,a,f16.6)") " Spatial extent of the function
:", sum(eigval) + + if (ifunc == 1) then + moment2nuc = 0 + do iatm = 1, ncenter + xtmp = a(iatm)%x - cenx + ytmp = a(iatm)%y - ceny + ztmp = a(iatm)%z - cenz + tmpval = a(iatm)%charge + moment2nuc(1, 1) = moment2nuc(1, 1) + xtmp*xtmp*tmpval + moment2nuc(2, 2) = moment2nuc(2, 2) + ytmp*ytmp*tmpval + moment2nuc(3, 3) = moment2nuc(3, 3) + ztmp*ztmp*tmpval + moment2nuc(1, 2) = moment2nuc(1, 2) + xtmp*ytmp*tmpval + moment2nuc(2, 3) = moment2nuc(2, 3) + ytmp*ztmp*tmpval + moment2nuc(1, 3) = moment2nuc(1, 3) + xtmp*ztmp*tmpval + end do + moment2nuc(3, 1) = moment2nuc(1, 3) + moment2nuc(2, 1) = moment2nuc(1, 2) + moment2nuc(3, 2) = moment2nuc(2, 3) + write (*, *) + write (*, "(' The quadrupole moment of nuclear charges:')") + write (*, "(' XX=',f16.8,' XY=',f16.8,' XZ=',f16.8)") moment2nuc(1, :) + write (*, "(' YX=',f16.8,' YY=',f16.8,' YZ=',f16.8)") moment2nuc(2, :) + write (*, "(' ZX=',f16.8,' ZY=',f16.8,' ZZ=',f16.8)") moment2nuc(3, :) + write (*, *) + write (*, "(' The quadrupole moment of the system:')") + write (*, "(' XX=',f16.8,' XY=',f16.8,' XZ=',f16.8)") moment2nuc(1, :) - moment2(1, :) + write (*, "(' YX=',f16.8,' YY=',f16.8,' YZ=',f16.8)") moment2nuc(2, :) - moment2(2, :) + write (*, "(' ZX=',f16.8,' ZY=',f16.8,' ZZ=',f16.8)") moment2nuc(3, :) - moment2(3, :) + end if + + else if (isel == 2 .or. isel == 5) then + realcenx = realcenx/intval + realceny = realceny/intval + realcenz = realcenz/intval + if (isel == 2) then + write (*, "(' Integral of the function:',1PE16.8,' a.u.')") intval + write (*, "(/,' Center of the function:')") + else + write (*, "(' Integral of the absolute of the function:',1PE16.8,' a.u.')") intval + write (*, "(/,' Center of the absolute of the function:')") + end if + write (*, "(' X=',f16.8,' Y=',f16.8,' Z=',f16.8,' Angstrom',/)") realcenx*b2a, realceny*b2a, realcenz*b2a + write (*, *) "Use this center for subsequent calculations? (y/n)" + read (*, *) selectyn + if (selectyn == 'y') then + cenx = realcenx + ceny = realceny + cenz = realcenz + end if + end if + end do +end subroutine -- 2.25.1