Fix gFortran compilation (gMultiWFN 3.7)

src/0123dim.f90

@@ -215,7 +215,7 @@ subroutine study1dim
 		    alive=.false.
 		    if (cubegenpath/=" ".and.ifiletype==1.and.isel==12) then
 			    inquire(file=cubegenpath,exist=alive)
-			    if (alive==.false.) then
+			    if (.not.alive) then
 				    write(*,"(a)") " Note: Albeit current file type is fch/fchk/chk and ""cubegenpath"" parameter in settings.ini has been defined, &
 				    the cubegen cannot be found, therefore electrostatic potential will still be calculated using internal code of Multiwfn"
 			    end if
@@ -2405,7 +2405,7 @@ subroutine setcontour
 		do while(.true.)
 			read(*,"(a)") extctrsetting
 			inquire(file=extctrsetting,exist=alive)
-			if (alive==.true.) exit
+			if (alive) exit
 			write(*,*) "File not found, input again"
 		end do
 		open(10,file=extctrsetting,status="old")

src/AdNDP.f90

@@ -667,7 +667,7 @@ subroutine AdNDP
 			write(*,*) "Note: If the suffix is .47, the Fock matrix will be directly loaded from it"
 			read(*,"(a)") c200tmp
 			inquire(file=c200tmp,exist=alive)
-			if (alive==.false.) then
+			if (.not.alive) then
 				write(*,*) "Error: Unable to find this file!"
 				cycle
 			end if

src/CDA.f90

@@ -1297,7 +1297,7 @@ subroutine plotintdiag(status,ifrag,jfrag,nCDAfrag,nmoCDA,FOcp,nmo0,nmo1,nmo2,oc
 		ilastplotted=0
 		cycle
 	end if
-	if ((idrawMObar==2.or.idrawMObar==3).and.FO1involveconn(iorb)==.false.) cycle
+	if ((idrawMObar==2.or.idrawMObar==3).and..not.FO1involveconn(iorb)) cycle
 	call solid !Use solid line to plot occupied orbital bars, use dashed line to plot virtual orbital bars
 	if (occ1(iorb)==0) call dash
 	call rline(xlow1,eneval,xhigh1,eneval)
@@ -1325,7 +1325,7 @@ subroutine plotintdiag(status,ifrag,jfrag,nCDAfrag,nmoCDA,FOcp,nmo0,nmo1,nmo2,oc
 		ilastplotted=0
 		cycle
 	end if
-	if ((idrawMObar==2.or.idrawMObar==3).and.FO2involveconn(iorb)==.false.) cycle
+	if ((idrawMObar==2.or.idrawMObar==3).and..not.FO2involveconn(iorb)) cycle
 	call solid
 	if (occ2(iorb)==0) call dash
 	call rline(xlow2,eneval,xhigh2,eneval)

src/DOS.f90

@@ -39,7 +39,7 @@ subroutine DOS
 character clegend*960 !Legend strings. (10+2) lines * 80 character per line
 integer :: legendx=400,legendy=160
 character :: TDOSstring*80="TDOS",OPDOSstring*80="OPDOS",graphformat_old*4
-character :: PDOSstring(nfragmax)*80=(/"PDOS frag.1","PDOS frag.2","PDOS frag.3","PDOS frag.4","PDOS frag.5","PDOS frag.6","PDOS frag.7","PDOS frag.8","PDOS frag.9","PDOS frag.10"/)
+character :: PDOSstring(nfragmax)*80=(/"PDOS frag.1 ","PDOS frag.2 ","PDOS frag.3 ","PDOS frag.4 ","PDOS frag.5 ","PDOS frag.6 ","PDOS frag.7 ","PDOS frag.8 ","PDOS frag.9 ","PDOS frag.10"/)
 integer :: ishowPDOSline(nfragmax),ishowPDOScurve(nfragmax),icurvewidth=3,ilinewidth=2,intarr(2)
 integer :: iclrPDOS(nfragmax)=(/ 1,3,10,14,12,9,13,11,6,7 /)
 !Below are used for defining fragments. For Mulliken/SCPA, they correspond to basis function, while for Hirshfeld/Becke, they correspond to atom indices
@@ -339,7 +339,7 @@ subroutine DOS
     read(*,"(a)") c200tmp
     if (c200tmp==" ") c200tmp="DOS.dat"
 	inquire(file=c200tmp,exist=alive)
-	if (alive==.false.) then
+	if (.not.alive) then
 	    write(*,*) "Error: Cannot find the file! Press ENTER button to return"
         read(*,*)
         cycle

src/GUI.f90

@@ -1037,7 +1037,7 @@ subroutine miniGUI
 !Show all orbitals
 subroutine showorbinfo(id)
 integer,intent (in) :: id
-character*3 :: orbtype(0:2)=(/ "A+B"," A"," B" /)
+character*3 :: orbtype(0:2)=(/ "A+B"," A "," B " /)
 character*6 :: symstr
 symstr=" "
 naorb=count(MOtype==1)
@@ -1066,7 +1066,7 @@ subroutine showorbinfo(id)
 !Show orbitals up to LUMO+10, works for wfntype==0,1,2
 subroutine showorbinfo2(id)
 integer,intent (in) :: id
-character*3 :: orbtype(0:2)=(/ "A+B"," A"," B" /)
+character*3 :: orbtype(0:2)=(/ "A+B"," A "," B " /)
 character*6 :: symstr
 symstr=" "
 naorb=count(MOtype==1)
@@ -1117,7 +1117,7 @@ subroutine showorbinfo2(id)
 !Show all occupied orbitals
 subroutine showorbinfo3(id)
 integer,intent (in) :: id
-character*3 :: orbtype(0:2)=(/ "A+B"," A"," B" /)
+character*3 :: orbtype(0:2)=(/ "A+B"," A "," B " /)
 character*6 :: symstr
 symstr=" "
 naorb=count(MOtype==1)
@@ -1461,7 +1461,7 @@ subroutine showorbsel(id,iorb)
 use function
 integer id,iorb
 real*8 molxlen,molylen,molzlen
-character*3 :: orbtype(0:2)=(/ "A+B"," A"," B" /)
+character*3 :: orbtype(0:2)=(/ "A+B"," A "," B " /)
 character*6 :: symstr
 ! Set grid for calculating cube data
 molxlen=(maxval(a%x)-minval(a%x))+2*aug3D

src/LSB.f90

@@ -846,7 +846,7 @@ subroutine integratebasinmix_LSB
 			rnowy=yarr(iy)
 			do ix=2,nx-1
 				rnowx=xarr(ix)
-				if (interbasgrid(ix,iy,iz)==.false.) cycle
+				if (.not.interbasgrid(ix,iy,iz)) cycle
  				nrefine=1
 				ndiv=nrefine**3
 				orgxref=rnowx-dx/2 !Take corner position as original point of microcycle

src/Makefile

@@ -1,11 +1,10 @@
-SIMD = -msse3
-OPT = -O2 -qopenmp -qopenmp-link=static -threads -qopt-matmul $(SIMD) -diag-disable 8290,8291,6371,10316 -fpp -mkl -static-intel
+OPT = -O2 -fopenmp -ffree-line-length-none -cpp
 #Options in the next line is for debugging purpose
 #OPT = -O0 -qopenmp -diag-disable 8290,8291,6371 -threads -qopenmp-link=static -debug all -g -traceback -check all -fstack-protector -fpp -mkl -static-intel
 
-LIB = ./dislin_d-11.0.a -lXm -lXt -lX11 -lGL #Works for CentOS 7.x
-FC = ifort
-EXE = Multiwfn
+LIB = ./dislin_d-11.0.a -lXm -lXt -lX11 -lGL -lblas -llapack #Works for CentOS 7.x
+FC = gfortran
+EXE = gMultiwfn
 LIBRETAPATH = ./libreta_slow
 
 objects = define.o util.o plot.o Bspline.o sym.o libreta.o function.o GUI.o sub.o integral.o Lebedev-Laikov.o \
@@ -168,7 +167,7 @@ ean.o: ${LIBRETAPATH}/ean.f90 hrr_012345.o eanvrr_012345.o boysfunc.o ${LIBRETAP
 	$(FC) $(OPT) -c ${LIBRETAPATH}/ean.f90
 
 hrr_012345.o: ${LIBRETAPATH}/hrr_012345.f90
-	$(FC) -O2 -diag-disable 6843 $(SIMD) -c ${LIBRETAPATH}/hrr_012345.f90
+	$(FC) $(OPT) -c ${LIBRETAPATH}/hrr_012345.f90
 
 eanvrr_012345.o: ${LIBRETAPATH}/eanvrr_012345.f90 boysfunc.o
 	$(FC) $(OPT) -c ${LIBRETAPATH}/eanvrr_012345.f90

src/Makefile_noGUI

@@ -1,9 +1,8 @@
-SIMD = -msse3
-OPT = -O2 -qopenmp -qopenmp-link=static -threads -qopt-matmul $(SIMD) -diag-disable 8290,8291,6371,10316 -fpp -mkl -static-intel
+OPT = -O2 -fopenmp -ffree-line-length-none -cpp
 
-LIB =
-FC = ifort
-EXE = Multiwfn
+LIB = -lblas -llapack
+FC = gfortran
+EXE = gMultiwfn
 LIBRETAPATH = ./libreta_slow
 
 objects = define.o util.o Bspline.o sym.o libreta.o function.o sub.o integral.o Lebedev-Laikov.o \
@@ -157,7 +156,7 @@ ean.o: ${LIBRETAPATH}/ean.f90 hrr_012345.o eanvrr_012345.o boysfunc.o ${LIBRETAP
 	$(FC) $(OPT) -c ${LIBRETAPATH}/ean.f90
 
 hrr_012345.o: ${LIBRETAPATH}/hrr_012345.f90
-	$(FC) -O2 -diag-disable 6843 $(SIMD) -c ${LIBRETAPATH}/hrr_012345.f90
+	$(FC) $(OPT) -c ${LIBRETAPATH}/hrr_012345.f90
 
 eanvrr_012345.o: ${LIBRETAPATH}/eanvrr_012345.f90 boysfunc.o
 	$(FC) $(OPT) -c ${LIBRETAPATH}/eanvrr_012345.f90

src/Multiwfn.f90

@@ -7,7 +7,7 @@ program multiwfn
 real*8,allocatable :: tmparr(:),tmparr2(:),tmpmat(:,:),tmpmat2(:,:) !For debug purpose
 integer,allocatable :: tmparri(:),tmparr2i(:),tmpmati(:,:),tmpmat2i(:,:)
 
-call kmp_set_warnings_off() !In rare case, "Cannot open message catalog "1041\libiomp5ui.dll"" may occurs, this calling avoid this problem, or user should set KMP_WARNINGS environment variable to 0
+!call kmp_set_warnings_off() !In rare case, "Cannot open message catalog "1041\libiomp5ui.dll"" may occurs, this calling avoid this problem, or user should set KMP_WARNINGS environment variable to 0
 
 !Try to get input file name from argument, which should be the first argument
 filename=" "
@@ -40,13 +40,13 @@ program multiwfn
 
 !For Linux/MacOS version, it seems the only way to set stacksize of each thread is to define KMP_STACKSIZE environment variable
 if (isys==1) then !Set via ompstacksize in settings.ini
-    call KMP_SET_STACKSIZE_S(ompstacksize)
+    !call KMP_SET_STACKSIZE_S(ompstacksize)
 else if (isys==2) then !The size should have been defined by KMP_STACKSIZE
     CALL getenv('KMP_STACKSIZE',c200tmp)
     if (c200tmp==" ") write(*,"(/,a)") " Warning: You should set ""KMP_STACKSIZE"" environment variable as mentioned in Section 2.1.2 of Multiwfn manual!"
 end if
 !write(*,"(' OpenMP stacksize for each thread: ',f10.2,' MB')") dfloat(KMP_GET_STACKSIZE_S())/1024/1024
-call mkl_set_num_threads(nthreads) !Use this to set number of cores used in MKL library (e.g. function matmul_blas)
+!call mkl_set_num_threads(nthreads) !Use this to set number of cores used in MKL library (e.g. function matmul_blas)
 
 write(*,*)
 
@@ -87,9 +87,9 @@ program multiwfn
 	end do
 	!Write current opened file to "lastfile" in settings.ini
 	inquire(file="settings.ini",exist=alive)
-	if (alive==.true.) then
+	if (alive) then
 		settingpath="settings.ini"
-	else if (alive==.false.) then
+	else if (.not.alive) then
 		call getenv("Multiwfnpath",c200tmp)
 		if (isys==1) then
 			settingpath=trim(c200tmp)//"\settings.ini"
@@ -543,7 +543,7 @@ program multiwfn
 		    else if (i==13) then
                 allocate(tmpmat(ncenter,ncenter))
                 inquire(file="bndmat.txt",exist=alive)
-	            if (alive==.false.) then
+	            if (.not.alive) then
 	                write(*,*) "Cannot find the bndmat.txt in current folder!"
                     cycle
                 end if

src/basin.f90

@@ -175,7 +175,7 @@ subroutine basinana
 			    if (igrid==1) then
                     continue
                 else if (igrid==2) then
-                    where (interbasgrid==0) cubmattmp=0
+                    where (.not.interbasgrid) cubmattmp=0
                 else if (igrid==3) then
                     if (ifuncbasin==1) then !The cubmat already records electron density
                         where(cubmat<0.001D0) cubmattmp=0
@@ -518,7 +518,7 @@ while the atomic indices correspond to the raw attractor indices. The atomic cha
 		write(*,"(' The number of unassigned grids:',i12)") numunassign
 		numgotobound=count(gridbas(2:nx-1,2:ny-1,2:nz-1)==-1)
 		write(*,"(' The number of grids travelled to box boundary:',i12)") numgotobound
-		where (grdposneg==.false.) cubmat=-cubmat !Recover original grid data
+		where (.not.grdposneg) cubmat=-cubmat !Recover original grid data
 		deallocate(grdposneg)
 
 		do iatt=1,numatt !Eliminate the attractors with very low value
@@ -551,7 +551,7 @@ while the atomic indices correspond to the raw attractor indices. The atomic cha
 		call clusdegenatt(0) !Cluster degenerate attractors as "real attractors" and calculate average coordinate and value for the real attractors
 
 		call detectinterbasgrd(6) !Detect interbasin grids
-		numinterbas=count(interbasgrid==.true.)
+		numinterbas=count(interbasgrid)
 		write(*,"(' The number of interbasin grids:',i12)") numinterbas
 
 	else if (isel==2) then
@@ -758,7 +758,7 @@ subroutine generatebasin(igridmethod)
 							attgrid(numatt,1)=inowx
 							attgrid(numatt,2)=inowy
 							attgrid(numatt,3)=inowz
-							if (grdposneg(inowx,inowy,inowz)==.true.) then
+							if (grdposneg(inowx,inowy,inowz)) then
 								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,cubmat(inowx,inowy,inowz)
 							else !This grid should has negative value
 								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,-cubmat(inowx,inowy,inowz)
@@ -859,7 +859,7 @@ subroutine generatebasin(igridmethod)
 							attgrid(numatt,1)=inowx
 							attgrid(numatt,2)=inowy
 							attgrid(numatt,3)=inowz
-							if (grdposneg(inowx,inowy,inowz)==.true.) then
+							if (grdposneg(inowx,inowy,inowz)) then
 								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,cubmat(inowx,inowy,inowz)
 							else !This grid should has negative value
 								write(*,"(i8,3f14.8,f20.8)") numatt,(orgx+(inowx-1)*dx)*b2a,(orgy+(inowy-1)*dy)*b2a,(orgz+(inowz-1)*dz)*b2a,-cubmat(inowx,inowy,inowz)
@@ -932,15 +932,15 @@ subroutine generatebasin(igridmethod)
 ! 		do itime=1,3 !Refine one time in general is sufficient
 		write(*,*) "Detecting boundary grids..."
 		call detectinterbasgrd(6) !It seems that using 26 directions to determine boundary grids doesn't bring evident benefit
-		write(*,"(' There are',i12,' grids at basin boundary')") count(interbasgrid==.true.)
+		write(*,"(' There are',i12,' grids at basin boundary')") count(interbasgrid)
 		write(*,*) "Refining basin boundary..."
 		!Below code is the adapted copy of above near-grid code
 		!$OMP PARALLEL DO private(ix,iy,iz,corrx,corry,corrz,ntrjgrid,inowx,inowy,inowz,valnow,imove,gradtmp,igradmax,gradmax,iatt,&
 		!$OMP ineiidx,idtmp,icorrx,icorry,icorrz,gradx,grady,gradz,sclgrad) shared(gridbas,attgrid) schedule(DYNAMIC) NUM_THREADS(nthreads)
 		do iz=2,nz-1
 			do iy=2,ny-1
 				do ix=2,nx-1
-					if (interbasgrid(ix,iy,iz)==.false.) cycle
+					if (.not.interbasgrid(ix,iy,iz)) cycle
 					if (gridbas(ix,iy,iz)<=0) cycle !Ignored the ones unassigned or gone to box boundary
 					ntrjgrid=0
 					inowx=ix
@@ -2067,7 +2067,7 @@ subroutine atmpopinbasin
 use basinintmod
 implicit real*8(a-h,o-z)
 inquire(file="basin.cub",exist=alive)
-if (alive==.false.) then
+if (.not.alive) then
 	write(*,*) "Error: basin.cub is not existed in current folder!"
 	return
 else
@@ -3329,7 +3329,7 @@ subroutine integratebasinmix(itype)
 		rnowy=yarr(iy)
 		do ix=2,nx-1
 			rnowx=xarr(ix)
-			if ((itype==2.or.itype==3).and.interbasgrid(ix,iy,iz)==.true.) cycle !If refine boundary grid at next stage, we don't calculate them at present stage
+			if ((itype==2.or.itype==3).and.interbasgrid(ix,iy,iz)) cycle !If refine boundary grid at next stage, we don't calculate them at present stage
 ! 			if (iz==nint(nz/2D0)) write(1,"('C',3f14.8)") rnowx*b2a,rnowy*b2a,rnowz*b2a !Examine grid distribution
 			iatt=gridbas(ix,iy,iz)
 ! 			do icp=1,numcp
@@ -3437,7 +3437,7 @@ subroutine integratebasinmix(itype)
 			rnowy=yarr(iy)
 			do ix=2,nx-1
 				rnowx=xarr(ix)
-				if (interbasgrid(ix,iy,iz)==.false.) cycle
+				if (.not.interbasgrid(ix,iy,iz)) cycle
 ! 				if (cubmat(ix,iy,iz)>0.001D0) then
 ! 					nrefine=2 !3 is the best
 ! 				else if (cubmat(ix,iy,iz)>0.001D0) then !0.0001 is the best
@@ -3597,7 +3597,7 @@ subroutine integratebasinmix(itype)
 	!$OMP END PARALLEL
 	call detectinterbasgrd(6)
 	write(*,*) "Basin boundary has been updated"
-	numinterbas=count(interbasgrid==.true.)
+	numinterbas=count(interbasgrid)
 	write(*,"(' The number of interbasin grids:',i12)") numinterbas	
 end if