1program multiwfn 2use defvar 3use util 4use function 5use topo 6implicit real*8(a-h,o-z) 7character nowdate*20,nowtime*20,inpstring*80,c200tmp*200,c200tmp2*200,c2000tmp*2000,outcubfile*200,selectyn,lovername*80,settingpath*200 8real*8 :: inx,iny,inz,tmpvec(3) 9integer :: iprintfunc=1 !The default function whose gradient and Hessian will be outputted at a point by main function 1 10integer,allocatable :: tmparrint(:) 11integer walltime1,walltime2 12real*8,allocatable :: d1add(:,:),d1min(:,:),d2add(:,:),d2min(:,:),d1addtmp(:,:),d1mintmp(:,:),d2addtmp(:,:),d2mintmp(:,:) !Store temporary data for drawing gradient map 13real*8,allocatable :: planemat_cust(:,:) !For storing temporary data of doing custom map 14real*8,allocatable :: planemat_bk(:,:) !Used to backup plane data 15real*8,allocatable :: tmpmat(:,:),tmparr(:) 16 17call getarg(1,filename) 18call getarg(2,cmdarg2) 1911 call loadsetting 20if (isys==1) write(*,*) "Multiwfn -- A Multifunctional Wavefunction Analyzer (for Windows 64bit)" 21if (isys==2) write(*,*) "Multiwfn -- A Multifunctional Wavefunction Analyzer (for Linux 64bit)" 22if (isys==3) write(*,*) "Multiwfn -- A Multifunctional Wavefunction Analyzer (for MacOS)" 23write(*,*) "Version 3.5(dev), release date: 2018-Jan-10" 24write(*,"(a)") " Project leader: Tian Lu (Beijing Kein Research Center for Natural Sciences)" 25write(*,*) "Citation of Multiwfn: Tian Lu, Feiwu Chen, J. Comput. Chem. 33, 580-592 (2012)" 26write(*,*) "Multiwfn official website: http://sobereva.com/multiwfn" 27write(*,*) "Multiwfn English forum: http://sobereva.com/wfnbbs" 28write(*,*) "Multiwfn Chinese forum: http://bbs.keinsci.com" 29 30!if (isys==1) call KMP_SET_STACKSIZE_S(ompstacksize) !For Linux/MacOS version, it seems the only way to set stacksize of each thread is to define KMP_STACKSIZE environment variable 31nthreads=getNThreads() 32call date_and_time(nowdate,nowtime) 33write(*,"(' ( The number of threads:',i3,' Current date: ',a,'-',a,'-',a,' Time: ',a,':',a,':',a,' )')") & 34nthreads,nowdate(1:4),nowdate(5:6),nowdate(7:8),nowtime(1:2),nowtime(3:4),nowtime(5:6) 35write(*,*) 36 37if (trim(filename)=="") then !Haven't defined filename variable 38 call mylover(lovername) 39 write(*,"(a,a,a)") " Input file path, for example E:\",trim(lovername),".wfn" 40 write(*,*) "(Supported: .wfn/.wfx/.fch/.molden/.31/.chg/.pdb/.xyz/.mol/.cub/.grd, etc.)" 41 write(*,"(a)") " Hint: Press ENTER button directly can select file in a GUI window. To reload the file last time used, simply input the letter ""o"". & 42 Input such as ?miku.fch can open the miku.fch in the same folder as the file last time used." 43 do while(.true.) 44 read(*,"(a)") filename 45 if (filename=='o') then 46 write(*,"(' The file last time used: ',a)") trim(lastfile) 47 filename=lastfile 48 end if 49 ltmp=len_trim(filename) 50 !Remove the first and the last " or 'symbol, because directly dragging file into the window will result in " or ' symbol, which is unrecognized by Multwifn 51 if (filename(1:1)=='"'.or.filename(1:1)=="'") filename(1:1)=" " 52 if (filename(ltmp:ltmp)=='"'.or.filename(ltmp:ltmp)=="'") filename(ltmp:ltmp)=" " 53 filename=adjustl(filename) 54 if (filename(1:1)=='?') then 55 do itmp=len_trim(lastfile),1,-1 56 if (isys==1.and.lastfile(itmp:itmp)=='\') exit 57 if ((isys==2.or.isys==3).and.lastfile(itmp:itmp)=='/') exit 58 end do 59 filename=lastfile(1:itmp)//trim(filename(2:)) 60 end if 61 inquire(file=filename,exist=alive) 62 if (alive.eqv..true.) exit 63 write(*,"('""',a,'"" ',a)") trim(filename),"cannot be found, input again" 64 end do 65 !Write current opened file to "lastfile" in settings.ini 66 inquire(file="settings.ini",exist=alive) 67 if (alive .eqv. .true.) then 68 settingpath="settings.ini" 69 else if (alive .eqv. .false.) then 70 call getenv("Multiwfnpath",c200tmp) 71 if (isys==1) then 72 settingpath=trim(c200tmp)//"\settings.ini" 73 else if (isys==2.or.isys==3) then 74 settingpath=trim(c200tmp)//"/settings.ini" 75 end if 76 end if 77 inquire(file=settingpath,exist=alive) 78 if (alive) then 79 open(20,file=settingpath,status="old") 80 call loclabel(20,"lastfile") 81 write(20,"(a)") "lastfile= "//trim(filename) 82 close(20) 83 end if 84else 85 inquire(file=filename,exist=alive) 86 if (alive.eqv..false.) then 87 write(*,*) "File not found, exit program..." 88 read(*,*) 89 stop 90 end if 91end if 92call readinfile(filename,0) 93write(*,"(/,3a)") " Loaded ",trim(filename)," successfully!" 94 95!!-- Backup various information of first loaded (meanwhile unmodified) molecule 96firstfilename=filename 97allocate(a_org(ncenter)) 98allocate(b_org(nprims)) 99allocate(CO_org(nmo,nprims)) 100allocate(MOocc_org(nmo)) 101a_org=a 102b_org=b 103CO_org=CO 104MOocc_org=MOocc 105nprims_org=nprims 106nmo_org=nmo 107ncenter_org=ncenter 108 109!!-- Initialize fragment 110nfragatmnum=ncenter !Default fragment is the whole molecule 111nfragatmnumbackup=ncenter 112allocate(fragatm(nfragatmnum),fragatmbackup(nfragatmnum)) 113forall (i=1:nfragatmnum) fragatm(i)=i 114forall (i=1:nfragatmnum) fragatmbackup(i)=i 115ifragcontri=0 116 117!!-- Call some routines only once 118if (ncenter>5000) write(*,"(a)") " Warning: There are very large number of many atoms, please wait very patiently for generating distance matrix..." 119call gendistmat !Generate distance matrix 120!Convert prebuild radii from Angstrom to Bohr. But some radii such as radii_hugo will remain unchanged since it is recorded as Bohr 121if (ifirstMultiwfn==1) then 122 vdwr=vdwr/b2a 123 vdwr_tianlu=vdwr_tianlu/b2a 124 covr=covr/b2a 125 covr_Suresh=covr_Suresh/b2a 126 covr_pyy=covr_pyy/b2a 127 covr_tianlu=covr_tianlu/b2a 128end if 129!Only get into dislin level 1 to get width and height of screen in pixels, don't do any other things 130!-- Show related molecular information 131if (ifiletype/=0.and.ifiletype/=8) then 132 call showformula 133 totmass=sum(atmwei(a%index)) 134 write(*,"(' Molecule weight:',f16.5)") totmass 135end if 136 137! call sys1eprop !Show some system 1e properties, only works when Cartesian basis functions are presented 138 139 140!!!--------------------- Now everything start ---------------------!!! 141!!!--------------------- Now everything start ---------------------!!! 142!!!--------------------- Now everything start ---------------------!!! 143do while(.true.) !Main loop 144 14510 write(*,*) 146if (allocated(cubmat)) write(*,*) "Note: A set of grid data presents in memory" 147write(*,*) " ------------ Main function menu ------------" 148! write(*,*) "-11 Load a new file" 149! write(*,*) "-10 Exit program" 150if (ifiletype/=7.and.ifiletype/=8) write(*,*) "0 Show molecular structure and view orbitals" 151if (ifiletype==7.or.ifiletype==8) write(*,*) "0 Show molecular structure and view isosurface" 152write(*,*) "1 Output all properties at a point" 153write(*,*) "2 Topology analysis" 154write(*,*) "3 Output and plot specific property in a line" 155write(*,*) "4 Output and plot specific property in a plane" 156write(*,*) "5 Output and plot specific property within a spatial region (calc. grid data)" 157write(*,*) "6 Check & modify wavefunction" 158write(*,*) "7 Population analysis and atomic charges" 159write(*,*) "8 Orbital composition analysis" 160write(*,*) "9 Bond order analysis" 161write(*,*) "10 Plot Total/Partial/Overlap population density-of-states (DOS)" 162write(*,*) "11 Plot IR/Raman/UV-Vis/ECD/VCD spectrum" 163write(*,*) "12 Quantitative analysis of molecular surface" 164if (allocated(cubmat)) write(*,*) "13 Process grid data" 165if (.not.allocated(cubmat)) write(*,*) "13 Process grid data (No grid data is presented currently)" 166write(*,*) "14 Adaptive natural density partitioning (AdNDP) analysis" 167write(*,*) "15 Fuzzy atomic space analysis" 168write(*,*) "16 Charge decomposition analysis (CDA) and extended CDA (ECDA)" 169write(*,*) "17 Basin analysis" 170write(*,*) "18 Electron excitation analysis" 171write(*,*) "100 Other functions (Part1)" 172write(*,*) "200 Other functions (Part2)" 173! write(*,*) "1000 Set some parameters" 174read(*,*) infuncsel1 175 176if (infuncsel1==-10) then !Exit program 177 stop 178else if (infuncsel1==-11) then !Load a new file 179 call dealloall 180 filename="" 181 deallocate(a_org,b_org,CO_org,MOocc_org,fragatm,fragatmbackup) 182 ifirstMultiwfn=0 183 goto 11 184end if 185 186!!! Every actual thing start from now on... 187 188!!!--------------------------------------- 189!1!!------------------- Show system structure and view isosurface of MOs or the grid data read from cube file 190if (infuncsel1==0) then ! 191 if (.not.(allocated(a).or.allocated(cubmat))) then 192 write(*,*) "Error: Data needed by this function is not presented! Check your input flie!" 193 write(*,*) "Press ENTER to continue" 194 read(*,*) 195 cycle 196 end if 197 if (ncenter>0) write(*,*) "Nucleus list:" 198 do i=1,ncenter 199 write(*,"(i5,'(',a2,')',' --> Charge:',f10.6,' x,y,z(Bohr):',3f11.6)") i,a(i)%name,a(i)%charge,a(i)%x,a(i)%y,a(i)%z 200 end do 201 if (allocated(CObasa).and.imodwfn==0) then !fch and occupation number hasn't been modified 202 if (wfntype==0) then 203 write(*,"(' Note: Orbital',i6,' is HOMO, energy:',f12.6,' a.u.',f12.6,' eV')") nint(nelec/2),MOene(nint(nelec/2)),MOene(nint(nelec/2))*au2eV 204 if (nint(nelec/2)+1<=nmo) then 205 write(*,"(' Orbital',i6,' is LUMO, energy:',f12.6' a.u.',f12.6,' eV')") nint(nelec/2)+1,MOene(nint(nelec/2)+1),MOene(nint(nelec/2)+1)*au2eV 206 gapene=MOene(nint(nelec/2)+1)-MOene(nint(nelec/2)) 207 write(*,"(' LUMO/HOMO gap:',f12.6,' a.u.',f12.6,' eV',f14.6,' kJ/mol')") gapene,gapene*au2eV,gapene*au2kJ 208 end if 209 else if (wfntype==1) then 210 write(*,"(' Range of alpha orbitals:',i5,' -',i5,' Range of Beta orbitals:',i5,' -',i5)") 1,nbasis,nbasis+1,nmo 211 write(*,"(' Note: Orbital',i6,' is HOMO of alpha spin, orbital',i6,' is HOMO of beta spin')") nint(naelec),nbasis+nint(nbelec) 212 if (nbasis>=nint(naelec)+1) then 213 gapenea=MOene(nint(naelec)+1)-MOene(nint(naelec)) 214 write(*,"(' LUMO/HOMO gap of alpha orbitals:',f12.6,' a.u.',f12.6,' eV')") gapenea,gapenea*au2eV 215 gapeneb=MOene(nbasis+nint(nbelec)+1)-MOene(nbasis+nint(nbelec)) 216 write(*,"(' LUMO/HOMO gap of beta orbitals: ',f12.6,' a.u.',f12.6,' eV')") gapeneb,gapeneb*au2eV 217 end if 218 else if (wfntype==2) then 219 write(*,"(' Index of SOMO orbitals:',7i6)") (i,i=nint(nbelec+1),nint(naelec)) 220 end if 221 end if 222 if (ifiletype==7.or.ifiletype==8) then !cube file 223 else 224 end if 225 226!!!--------------------------------------- 227!1!!------------------- Output properties at a point 228else if (infuncsel1==1) then 229 do while(.true.) 230 write(*,"(a)") " Input x,y,z, divided by space or comma, e.g. 3.3,2.0,-0.3" 231 write(*,"(a)") " or input e.g. ""a5"" to use nuclear position of atom 5" 232 write(*,"(a)") " input e.g. ""o8"" to select orbital 8, whose wavefunction value will be shown" 233 write(*,"(a)") " input e.g. ""f3"" to select function 3, whose gradient and Hessian will be shown, input ""allf"" can print all available functions" 234 write(*,"(a)") " input ""d"" to decompose a property at a point to orbital contributions" 235 write(*,"(a)") " input ""q"" to return" 236 read(*,"(a)") inpstring 237 inpstring=adjustl(inpstring) 238 if (inpstring(1:1)=='q') then 239 exit 240 else if (inpstring(1:4)=='allf') then 241 write(*,*) "1 Electron density (Analytical Hessian)" 242 write(*,*) "3 Laplacian of electron density" 243 write(*,*) "4 Value of orbital wavefunction" 244 if (ELFLOL_type==0) write(*,*) "9 Electron localization function(ELF)" 245 if (ELFLOL_type==1) write(*,*) "9 Electron localization function(ELF) defined by Tsirelson" 246 if (ELFLOL_type==2) write(*,*) "9 Electron localization function(ELF) defined by Lu, Tian" 247 if (ELFLOL_type==0) write(*,*) "10 Localized orbital locator(LOL)" 248 if (ELFLOL_type==1) write(*,*) "10 Localized orbital locator(LOL) defined by Tsirelson" 249 if (ELFLOL_type==2) write(*,*) "10 Localized orbital locator(LOL) defined by Lu, Tian" 250 write(*,*) "12 Total electrostatic potential" 251 write(*,*) "100 User-defined real space function" 252 else if (inpstring(1:1)=='f') then 253 read(inpstring(2:),*) iprintfunc 254 write(*,"(' Real space function',i5,' is selected')") iprintfunc 255 else if (inpstring(1:1)=='a') then 256 read(inpstring(2:),*) iatm 257 if (iatm>0.and.iatm<=ncenter) then 258 write(*,*) "Note: Unless otherwise specified, all units are in a.u." 259 write(*,"(' Atom',i6,' X,Y,Z(Bohr):',3f14.8)") iatm,a(iatm)%x,a(iatm)%y,a(iatm)%z 260 call showptprop(a(iatm)%x,a(iatm)%y,a(iatm)%z,iprintfunc,6) 261 else 262 write(*,*) "The index exceeds valid range" 263 end if 264 else if (inpstring(1:1)=='o') then 265 read(inpstring(2:),*) iorbseltmp 266 if (iorbseltmp>0.and.iorbseltmp<=nmo) then 267 iorbsel=iorbseltmp 268 else 269 write(*,*) "The index exceeds valid range" 270 end if 271 else if (inpstring(1:1)=='d') then 272 write(*,*) "Input x,y,z of the point, e.g. 3.3,2.0,-0.3" 273 read(*,*) inx,iny,inz 274 write(*,*) "You inputted coordinate is in which unit? 1:Bohr 2:Angstrom" 275 read(*,*) iunit 276 if (iunit==2) then 277 inx=inx/b2a 278 iny=iny/b2a 279 inz=inz/b2a 280 end if 281 call decompptprop(inx,iny,inz) 282 else 283 read(inpstring,*) inx,iny,inz 284 write(*,*) "You inputted coordinate is in which unit? 1:Bohr 2:Angstrom" 285 read(*,*) iunit 286 if (iunit==2) then 287 inx=inx/b2a 288 iny=iny/b2a 289 inz=inz/b2a 290 end if 291 write(*,*) "Note: Unless otherwise specified, all units are in a.u." 292 call showptprop(inx,iny,inz,iprintfunc,6) 293 end if 294 write(*,*) 295 end do 296 297!!!--------------------------------------- 298!2!!------------------- Topology analysis 299else if (infuncsel1==2) then 300 call delvirorb(1) !Don't need virtual orbitals, delete them for faster calculation 301 call topomain 302 303 304 305 306!!!--------------------------------------- 307!!!--------------------------------------- 308!3!------------------- Draw property in a line 309!!!--------------------------------------- 310!!!--------------------------------------- 311else if (infuncsel1==3) then 312 ncustommap=0 !Clean custom operation setting that possibly defined by other modules 313 if (allocated(custommapname)) deallocate(custommapname) 314 if (allocated(customop)) deallocate(customop) 315 write(*,*) "-10 Return to main menu" 316 write(*,*) "-2 Obtain deformation property" 317 write(*,*) "-1 Obtain promolecule property" 318 write(*,*) "0 Set custom operation" 319300 call selfunc_interface(infuncsel2) 320 321 if (infuncsel2==0.or.infuncsel2==-1.or.infuncsel2==-2) then 322 if (infuncsel2==0) call customplotsetup 323 if (infuncsel2==-1) then 324 ipromol=1 325 else 326 ipromol=0 327 end if 328 if (infuncsel2==-1.or.infuncsel2==-2) call setPromol 329 write(*,*) "-10 Return to main menu" 330 goto 300 331 else if (infuncsel2==-10) then 332 cycle 333 else if (infuncsel2==111) then 334 write(*,*) "Input indices of two atoms, e.g. 1,4, or input an atom and zero, e.g. 5,0" 335 read(*,*) iatmbecke1,iatmbecke2 336 end if 337 338301 write(*,"(a,f8.4,a)") " 0 Set extension distance for mode 1, current:",aug1D," Bohr" 339 write(*,*) "1 Input index of two nuclei to define a line" 340 write(*,*) "2 Input coordinate of two points to define a line" 341 read(*,*) infuncsel3 342 343 if (infuncsel3==0) then 344 write(*,*) "Input augment distance (in Bohr, e.g. 2.5)" 345 read(*,*) aug1D 346 goto 301 347 else if (infuncsel3==1) then 348 do while(.true.) 349 write(*,*) "Input two number to select two nuclei (e.g. 1,3)" 350 read(*,*) iselatm1,iselatm2 351 if (iselatm1/=iselatm2.and.min(iselatm1,iselatm2)>=1.and.max(iselatm1,iselatm2)<=ncenter) exit 352 write(*,*) "Invalid input" 353 end do 354 write(*,"(' Nucleus',i5,' Charge:',f6.2,' X,Y,Z:',3f12.6)") iselatm1,a(iselatm1)%charge,a(iselatm1)%x,a(iselatm1)%y,a(iselatm1)%z 355 write(*,"(' Nucleus',i5,' Charge:',f6.2,' X,Y,Z:',3f12.6)") iselatm2,a(iselatm2)%charge,a(iselatm2)%x,a(iselatm2)%y,a(iselatm2)%z 356 torgx=a(iselatm1)%x 357 torgy=a(iselatm1)%y 358 torgz=a(iselatm1)%z 359 tendx=a(iselatm2)%x 360 tendy=a(iselatm2)%y 361 tendz=a(iselatm2)%z 362 ratio=dsqrt((tendx-torgx)**2+(tendy-torgy)**2+(tendz-torgz)**2)/aug1D 363 orgx1D=torgx-(tendx-torgx)/ratio 364 orgy1D=torgy-(tendy-torgy)/ratio 365 orgz1D=torgz-(tendz-torgz)/ratio 366 endx1D=tendx+(tendx-torgx)/ratio 367 endy1D=tendy+(tendy-torgy)/ratio 368 endz1D=tendz+(tendz-torgz)/ratio 369 totdist=dsqrt((endx1D-orgx1D)**2+(endy1D-orgy1D)**2+(endz1D-orgz1D)**2) 370 atomr1=aug1D 371 atomr2=totdist-aug1D 372 else if (infuncsel3==2) then 373 write(*,*) "Input x1,y1,z1,x2,y2,z2 to define two points (in Bohr)" 374 write(*,*) "e.g. 0,0,3.2,-1,-0.26,2.8" 375 read(*,*) x1,y1,z1,x2,y2,z2 376 orgx1D=x1 377 orgy1D=y1 378 orgz1D=z1 379 endx1D=x2 380 endy1D=y2 381 endz1D=z2 382 atomr1=0D0 383 atomr2=0D0 384 end if 385 npointcurve=num1Dpoints !The number of data to plot 386 if (infuncsel2==12) then 387 npointcurve=num1Dpoints/6 !Calculate ESP is time consuming, so decrease the number of points 388 write(*,*) "Please wait..." 389 end if 390 if (infuncsel2/=12.and.expcutoff<0) write(*,"(' Note: All exponential functions exp(x) with x<',f8.3,' is be ignored ')") expcutoff 391 392 if (allocated(curvex)) deallocate(curvex) 393 if (allocated(curvey)) deallocate(curvey) 394 allocate(curvex(npointcurve)) 395 allocate(curvey(npointcurve)) 396 if (allocated(curveytmp)) deallocate(curveytmp) 397 if (ncustommap/=0) allocate(curveytmp(npointcurve)) 398 transx=(endx1D-orgx1D)/npointcurve 399 transy=(endy1D-orgy1D)/npointcurve 400 transz=(endz1D-orgz1D)/npointcurve 401 transr=dsqrt(transx**2+transy**2+transz**2) 402 403 write(*,*) 404 write(*,"(' Original point in X,Y,Z: ',3f10.5)") orgx1D,orgy1D,orgz1D !Output grid information 405 write(*,"(' End point in X,Y,Z: ',3f10.5)") endx1D,endy1D,endz1D 406 write(*,"(' Translation vector in X,Y,Z:',3f10.5,' Norm:',f10.5)") transx,transy,transz,transr 407 write(*,"(' Number of points:',i10)") npointcurve 408 409 icustom=0 410 curvey=0D0 411 if (ipromol==1) goto 311 412310 continue 413nthreads=getNThreads() 414!$OMP parallel do shared(curvex,curvey) private(i,rnowx,rnowy,rnowz) num_threads(nthreads) 415 do i=1,npointcurve !Calculate data for line plot 416 rnowx=orgx1D+i*transx 417 rnowy=orgy1D+i*transy 418 rnowz=orgz1D+i*transz 419 curvex(i)=i*transr 420 if (infuncsel2==111) then 421 curvey(i)=beckewei(rnowx,rnowy,rnowz,iatmbecke1,iatmbecke2) 422 else 423 curvey(i)=calcfuncall(infuncsel2,rnowx,rnowy,rnowz) 424 end if 425 end do 426!$OMP end parallel do 427 428311 if (ncustommap/=0) then !Calculate data for custom map 429 if (icustom==0) then 430 curveytmp=curvey 431 else if (icustom/=0) then 432 if (customop(icustom)=='+') curveytmp=curveytmp+curvey 433 if (customop(icustom)=='-') curveytmp=curveytmp-curvey 434 if (customop(icustom)=='x'.or.customop(icustom)=='*') curveytmp=curveytmp*curvey 435 if (customop(icustom)=='/') curveytmp=curveytmp/curvey 436 end if 437 if (icustom/=ncustommap) then 438 icustom=icustom+1 439 filename=custommapname(icustom) 440 call dealloall 441 write(*,"(' Loading: ',a)") trim(filename) 442 call readinfile(filename,1) 443 !Generate temporary fragatm 444 deallocate(fragatm) 445 nfragatmnum=ncenter 446 allocate(fragatm(nfragatmnum)) 447 do iatm=1,ncenter 448 fragatm(iatm)=iatm 449 end do 450 !Input the MO index for current file. Since the MO index may be not the same as the first loaded one 451 if (infuncsel2==4) then 452 write(*,"(' Input the index of the orbital to be calculated for ',a,' e.g. 3')") trim(filename) 453 read(*,*) iorbsel 454 end if 455 goto 310 456 else 457 curvey=curveytmp 458 call dealloall 459 write(*,"(' Reloading: ',a)") trim(firstfilename) 460 call readinfile(firstfilename,1) 461 !Recovery user-defined fragatm from the backup 462 deallocate(fragatm) 463 nfragatmnum=nfragatmnumbackup 464 allocate(fragatm(nfragatmnum)) 465 fragatm=fragatmbackup 466 end if 467 end if 468 469 write(*,"(' Minimal/Maximum value:',2D16.8)") minval(curvey),maxval(curvey) 470 write(*,"(' Summing up all values:',D18.8,' Integration value:',D18.8)") sum(curvey),sum(curvey)*transr 471 exty=(maxval(curvey)-minval(curvey))/10 472 if (exty<maxval(abs(curvey))/100) exty=maxval(abs(curvey))/10 !Sometimes the difference between maximal and minimal values are too small, so do special treatment 473 curveymin=minval(curvey)-exty 474 curveymax=maxval(curvey)+exty 475 steplabx=maxval(curvex)/10 476 steplaby=exty 477 i=-1 478 479 do while(.true.) 480 if (isilent==0.and.i==-1) then 481 if (atomr1==atomr2) then 482 else !Draw the two atom positions 483 end if 484 end if 485 write(*,*) 486 write(*,*) "-1 Show the graph again" 487 write(*,*) "0 Return to main menu" 488 write(*,*) "1 Save the graph to a file in current folder" 489 write(*,*) "2 Export the data to line.txt in current folder" 490 write(*,"(a,1PE14.6,a,1PE14.6)") " 3 Change range of Y axis, current: from",curveymin," to",curveymax 491 if (icurve_vertlinex==0) write(*,*) "4 Draw a vertical line with specific X" 492 if (icurve_vertlinex==1) write(*,*) "4 Delete the vertical line" 493 write(*,"(' 5 Change the ratio of X and Y length, current:',f10.5)") curvexyratio 494 write(*,"(' 6 Find the positions of local minimum and maximum')") 495 write(*,"(' 7 Find the positions where function value equals to specified value')") 496 if (ilog10y==0) write(*,*) "8 Use logarithmic scaling of Y axis" 497 if (ilog10y==1) write(*,*) "8 Use linear scaling of Y axis" 498 write(*,*) "9 Change the line color" 499 if (ilog10y==0) write(*,"(a,f8.3,1PE14.5)") " 10 Set stepsize in X and Y axes, current:",steplabx,steplaby 500 if (ilog10y==1) write(*,"(a,f8.3)") " 10 Set the stepsize in X axis, current:",steplabx 501 if (ilenunit1D==1) write(*,*) "11 Change length unit of the graph to Angstrom" 502 if (ilenunit1D==2) write(*,*) "11 Change length unit of the graph to Bohr" 503 write(*,"(a,i3)") " 12 Set width of curve line, current:",icurvethick 504 505 read(*,*) i 506 if (i==-1) then 507 cycle 508 else if (i==0) then 509 exit 510 else if (i==1) then 511 write(*,"(a,a,a)") " Graph have been saved to ",trim(graphformat)," file with ""DISLIN"" prefix in current directory" 512 else if (i==2) then 513 open(14,file="line.txt",status="replace") 514 do i=1,npointcurve !Output result to file 515 rnowx=orgx1D+i*transx 516 rnowy=orgy1D+i*transy 517 rnowz=orgz1D+i*transz 518 curvex(i)=i*transr 519 write(14,"(4f12.6,1PE18.10)") rnowx*b2a,rnowy*b2a,rnowz*b2a,curvex(i)*b2a,curvey(i) 520 end do 521 close (14) 522 write(*,*) "Results have been saved to line.txt in current folder" 523 write(*,"(a)") " Unit is Angstrom. The first three columns are actual coordinates, the fourth column & 524 is X position in the curve graph, the fifth column is function value" 525 else if (i==3) then 526 if (ilog10y==0) write(*,*) "Input minimum and maximum value of Y axis e.g. -0.1,2" 527 if (ilog10y==1) write(*,*) "Input minimum and maximum value of Y axis e.g. -1,5 means from 10^-1 to 10^5" 528 read(*,*) curveymin,curveymax 529 else if (i==4.and.icurve_vertlinex==0) then 530 write(*,*) "Input X (in Bohr) e.g. 1.78" 531 read(*,*) curve_vertlinex 532 icurve_vertlinex=1 533 else if (i==4.and.icurve_vertlinex==1) then 534 icurve_vertlinex=0 535 else if (i==5) then 536 write(*,*) "Input a value" 537 read(*,*) curvexyratio 538 else if (i==6) then 539 numlocmin=0 540 numlocmax=0 541 do ipoint=2,npointcurve-1 542 gradold=curvey(ipoint)-curvey(ipoint-1) 543 gradnew=curvey(ipoint+1)-curvey(ipoint) 544 if (gradold*gradnew<0D0) then 545 if (gradold>gradnew) then 546 numlocmax=numlocmax+1 547 write(*,"(' Local maximum X (Bohr):',f12.6,' Value:',D18.8)") curvex(ipoint),curvey(ipoint) 548 else if (gradold<gradnew) then 549 numlocmin=numlocmin+1 550 write(*,"(' Local minimum X (Bohr):',f12.6,' Value:',D18.8)") curvex(ipoint),curvey(ipoint) 551 end if 552 end if 553 end do 554 write(*,"(' Totally found',i5,' local minimum,',i5,' local maximum')") numlocmin,numlocmax 555 else if (i==7) then 556 write(*,*) "Input a value" 557 read(*,*) specvalue 558 numfind=0 559 do ipoint=1,npointcurve-1 560 if ( (specvalue>curvey(ipoint).and.specvalue<curvey(ipoint+1)) .or.& 561 (specvalue<curvey(ipoint).and.specvalue>curvey(ipoint+1)) ) then 562 !Use linear interpolation to evaluate the X position 563 tmpratio=abs(specvalue-curvey(ipoint))/abs(curvey(ipoint+1)-curvey(ipoint)) 564 specvaluex=curvex(ipoint)+tmpratio*transr 565 numfind=numfind+1 566 write(*,"(' #',i5,' X (Bohr):',f12.6)") numfind,specvaluex 567 end if 568 end do 569 if (numfind==0) write(*,*) "Found nothing" 570 else if (i==8) then 571 if (ilog10y==1) then 572 ilog10y=0 573 !Recover default limit 574 curveymin=minval(curvey)-(maxval(curvey)-minval(curvey))/10 575 curveymax=maxval(curvey)+(maxval(curvey)-minval(curvey))/10 576 else if (ilog10y==0) then 577 ilog10y=1 578 write(*,*) "Input minimum and maximum value of Y axis e.g. -1,5 means from 10^-1 to 10^5" 579 read(*,*) curveymin,curveymax 580 end if 581 else if (i==9) then 582 write(*,*) "Use which color?" 583 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 584 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 585 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 586 read(*,*) iclrcurve 587 else if (i==10) then 588 if (ilog10y==0) then 589 write(*,*) "Input the step size between the labels in X and Y axes, respectively" 590 write(*,*) "e.g. 1.5,20" 591 read(*,*) steplabx,steplaby 592 else if (ilog10y==1) then 593 write(*,*) "Input the step size between the labels in X axis, e.g. 1.5" 594 read(*,*) steplabx 595 end if 596 else if (i==11) then 597 if (ilenunit1D==1) then 598 ilenunit1D=2 !Angstrom 599 else if (ilenunit1D==2) then 600 ilenunit1D=1 !Bohr 601 end if 602 else if (i==12) then 603 write(*,*) "Input line width, e.g. 5" 604 read(*,*) icurvethick 605 end if 606 end do 607 608 609 610!!!---------------------------------------- 611!!!---------------------------------------- 612!4!----------------------- Draw plane graph 613!!!---------------------------------------- 614!!!---------------------------------------- 615else if (infuncsel1==4) then 616 ncustommap=0 !Clean custom operation setting that possibly defined by other modules 617 if (allocated(custommapname)) deallocate(custommapname) 618 if (allocated(customop)) deallocate(customop) 619 if (allocated(tmparrint)) deallocate(tmparrint) 620 write(*,*) "-10 Return to main menu" 621 write(*,*) "-2 Obtain of deformation property" 622 write(*,*) "-1 Obtain of promolecule property" 623 write(*,*) "0 Set custom operation" 624400 call funclist 625 read(*,*) infuncsel2 626 if (infuncsel2==4) then !4=MO value, 111=becke weighting function 627 write(*,"(a,i10)") " Input the orbital index that you want to plot, should between 1 and",nmo 628 write(*,"(a)") " If you want to plot contour map for two orbitals at the same time, input two indices and separate them by comma, e.g. 8,10" 629 read(*,"(a)") c200tmp 630 if (index(c200tmp,',')/=0) then !Inputted two orbitals 631 read(c200tmp,*) iorbsel,iorbsel2 632 else 633 read(c200tmp,*) iorbsel 634 end if 635 else if (infuncsel2==20) then !Input length scale to evaluate EDR(r;d) 636 write(*,*) "The EDR(r;d) computing code was contributed by Arshad Mehmood" 637 write(*,"(a,/)") " References: J. Chem. Phys., 141, 144104 (2014); J. Chem. Theory Comput., 12, 79 (2016); Angew. Chem. Int. Ed., 56, 6878 (2017)" 638 write(*,*) " Input length scale d (Bohr) e.g. 0.85" 639 read(*,*) dedr 640 else if (infuncsel2==21) then !Input parameters to evaluate D(r) 641 write(*,*) "The D(r) computing code was contributed by Arshad Mehmood" 642 write(*,"(a,/)") " References: J. Chem. Theory Comput., 12, 3185 (2016); Phys. Chem. Chem. Phys., 17, 18305 (2015)" 643 write(*,*) "1 Manually input total number, start and increment in EDR exponents" 644 write(*,*) "2 Use default values i.e. 20,2.50,1.50" 645 read(*,*) edrmaxpara 646 if (edrmaxpara==1) then 647 write(*,*) "Please input in order: exponents start increment e.g. 20,2.5,1.5" 648 write(*,*) "Note: Max. allowed exponents are 50 and min. allowed increment is 1.01" 649 read(*,*) nedr,edrastart,edrainc 650 if (nedr<1) then 651 write(*,*) "Error: Bad Number of EDR exponents. Should be between 1 to 50" 652 write(*,*) "Press ENTER to exit" 653 read(*,*) 654 stop 655 else if (nedr>50) then 656 write(*,*) "Error: Bad Number of EDR exponents. Should be between 1 to 50" 657 write(*,*) "Press ENTER to exit" 658 read(*,*) 659 stop 660 end if 661 if (edrainc<1.01d0) then 662 write(*,*) "Error: Bad increment in EDR exponents. Should not be less than 1.01" 663 write(*,*) "Press ENTER to exit" 664 read(*,*) 665 stop 666 end if 667 else if (edrmaxpara==2) then 668 nedr=20 669 edrastart=2.5d0 670 edrainc=1.5d0 671 end if 672 write(*,*) "The following EDR exponents will be used in calculation:" 673 wrtstart=edrastart 674 do wrtnumedr=1,nedr 675 wrtexpo(wrtnumedr)=wrtstart 676 wrtstart=wrtstart/edrainc 677 write(*,"(E13.5)") wrtexpo(wrtnumedr) 678 end do 679 write(*,*) 680 else if (infuncsel2==111) then !Calculate Becke weighting function 681 write(*,*) "Input indices of two atoms to calculate Becke overlap weight, e.g. 1,4" 682 write(*,*) "or input index of an atom and zero to calculate Becke atomic weight, e.g. 5,0" 683 read(*,*) iatmbecke1,iatmbecke2 684 else if (infuncsel2==112) then !Calculate Hirshfeld weighting function 685 write(*,*) "Input index of the atoms you are interested in, e.g. 2,3,7-10" 686 read(*,"(a)") c2000tmp 687 call str2arr(c2000tmp,ntmp) 688 allocate(tmparrint(ntmp)) 689 call str2arr(c2000tmp,ntmp,tmparrint) 690 write(*,"(a)") " How to generate the atomic densities that used in the calculation of Hirshfeld weight?" 691 write(*,*) "1 Based on atomic .wfn files" 692 write(*,*) "2 Based on built-in atomic densities (see Appendix 3 of the manual for detail)" 693 read(*,*) iHirshdenstype 694 if (iHirshdenstype==1) call setpromol 695 else if (infuncsel2==500.or.infuncsel2==510.or.infuncsel2==511.or.infuncsel2==512) then !Calculate rho(A)*ln[rho(A)/rho0(A)], or rho(A), or rho0(A) 696 call setpromol 697 allocate(tmparrint(1)) 698 write(*,*) "Input index of the atom you are interested in, e.g. 4" 699 read(*,*) tmparrint 700 else if (infuncsel2==501.or.infuncsel2==502.or.infuncsel2==503) then !sum[rhoA*ln(rhoA/rho0A), sum[(rhoA-rho0A)/rhoA], difference between relative entropy and deformation density 701 call setpromol 702 end if 703 704 inucespplot=0 705 if (infuncsel2==8) inucespplot=1 !For deal with plotting nucesp property, special treatment is needed 706 imarkrefpos=0 707 if (infuncsel2==17.or.infuncsel2==19.or.(infuncsel2==100.and.iuserfunc==24)) imarkrefpos=1 !Only for correlation hole/factor, source function, linear response kernel... reference point will be marked on contour map 708 709 if (infuncsel2==0.or.infuncsel2==-1.or.infuncsel2==-2) then 710 if (infuncsel2==0) call customplotsetup 711 if (infuncsel2==-1) then 712 ipromol=1 713 else 714 ipromol=0 715 end if 716 if (infuncsel2==-1.or.infuncsel2==-2) call setPromol 717 write(*,*) "-10 Return to main menu" 718 goto 400 719 else if (infuncsel2==-10) then 720 cycle 721 end if 722 723 if (iorbsel2/=0) then 724 idrawtype=2 !Only contour line map is available when another orbital is needed to be plotted together 725 else if (iorbsel2==0) then 726 write(*,*) " -10 Return to main menu" 727 write(*,*) "Draw which type of map?" 728 write(*,*) "1 Color-filled map" 729 write(*,*) "2 Contour line map" 730 write(*,*) "3 Relief map" 731 write(*,*) "4 Shaded surface map" 732 write(*,*) "5 Shaded surface map with projection" 733 write(*,*) "6 Gradient lines map with/without contour lines" 734 write(*,*) "7 Vector field map with/without contour lines" 735 read(*,*) idrawtype 736 end if 737 738 if (idrawtype==1.or.idrawtype==2.or.idrawtype==6.or.idrawtype==7) then !Initialize contour line setting 739 if (idrawtype==2.or.idrawtype==6) idrawcontour=1 740 if (idrawtype==1.or.idrawtype==7) idrawcontour=0 741 if (idrawtype/=1) iatom_on_contour=1 742 ilabel_on_contour=0 743 ctrval=0.0D0 744 if (infuncsel2==9.or.infuncsel2==10) then !A special contour setting suitable for ELF and LOL 745 lastctrval=21 746 do i=1,21 747 ctrval(i)=(i-1)*0.05D0 748 end do 749 else !General contour setting for other real space functions 750 lastctrval=62 751 ctrval(1)=1D-3 752 do i=0,9 !set the value of contour line 753 ctrval(3*i+2)=2*1D-3*10**i 754 ctrval(3*i+3)=4*1D-3*10**i 755 ctrval(3*i+4)=8*1D-3*10**i 756 end do 757 ctrval(32:62)=-ctrval(1:31) 758 end if 759 else if (idrawtype==-10) then 760 cycle 761 end if 762 763 write(*,*) " -10 Return to main menu" 764 write(*,*) "How many grids in the two dimensions respectively?" 765 if (idrawtype==1.or.idrawtype==2.or.idrawtype==6) then 766 if (infuncsel2==12.or.(infuncsel2==100.and.iuserfunc==60).or.iuserfunc==39.or.iuserfunc==101.or.iuserfunc==102) then 767 write(*,*) "(100,100 is recommended)" !Because calculating ESP is very time consuming, so use lower grid 768 else 769 write(*,*) "(200,200 is recommended)" 770 end if 771 else if (idrawtype==3.or.idrawtype==4.or.idrawtype==5) then 772 write(*,*) "(100,100 is recommended)" 773 else if (idrawtype==7) then 774 write(*,*) "(80,80 is recommended)" 775 end if 776 write(*,*) "Hint: You can press ENTER button directly to use recommended value" 777 read(*,"(a)") c200tmp 778 if (c200tmp==' ') then !Press enter directly 779 if (idrawtype==1.or.idrawtype==2.or.idrawtype==6) then 780 if (infuncsel2==12.or.(infuncsel2==100.and.iuserfunc==60).or.iuserfunc==39.or.iuserfunc==101.or.iuserfunc==102) then 781 ngridnum1=100 782 ngridnum2=100 783 else 784 ngridnum1=200 785 ngridnum2=200 786 end if 787 else if (idrawtype==3.or.idrawtype==4.or.idrawtype==5) then 788 ngridnum1=100 789 ngridnum2=100 790 else if (idrawtype==7) then 791 ngridnum1=80 792 ngridnum2=80 793 end if 794 else if (c200tmp=='-10') then 795 cycle 796 else 797 read(c200tmp,*) ngridnum1,ngridnum2 798 end if 799 800 if (allocated(planemat)) deallocate(planemat) 801 if (allocated(planemattmp)) deallocate(planemattmp) 802 allocate(planemat(ngridnum1,ngridnum2),planemattmp(ngridnum1,ngridnum2)) !planemattmp is used in many cases below 803 if (ncustommap/=0) then 804 if (allocated(planemat_cust)) deallocate(planemat_cust) 805 if (allocated(d1addtmp)) deallocate(d1addtmp) 806 if (allocated(d1mintmp)) deallocate(d1mintmp) 807 if (allocated(d2addtmp)) deallocate(d2addtmp) 808 if (allocated(d2mintmp)) deallocate(d2mintmp) 809 allocate(planemat_cust(ngridnum1,ngridnum2)) 810 allocate(d1addtmp(ngridnum1,ngridnum2)) 811 allocate(d1mintmp(ngridnum1,ngridnum2)) 812 allocate(d2addtmp(ngridnum1,ngridnum2)) 813 allocate(d2mintmp(ngridnum1,ngridnum2)) 814 end if 815 if (idrawtype==6.or.idrawtype==7) then !Draw gradient lines 816 if (allocated(d1add)) deallocate(d1add) 817 if (allocated(d1min)) deallocate(d1min) 818 if (allocated(d2add)) deallocate(d2add) 819 if (allocated(d2min)) deallocate(d2min) 820 if (allocated(gradd1)) deallocate(gradd1) 821 if (allocated(gradd2)) deallocate(gradd2) 822 allocate(d1add(ngridnum1,ngridnum2)) 823 allocate(d1min(ngridnum1,ngridnum2)) 824 allocate(d2add(ngridnum1,ngridnum2)) 825 allocate(d2min(ngridnum1,ngridnum2)) 826 allocate(gradd1(ngridnum1,ngridnum2)) 827 allocate(gradd2(ngridnum1,ngridnum2)) 828 end if 829 830 write(*,*) " -10 Return to main menu" 831401 write(*,*) "Define the plane to be plotted" 832 write(*,*) "1:XY 2:XZ 3:YZ 4:Define by three atoms 5:Define by three points" 833 write(*,*) "6:Input origin and translation vector (For expert user)" 834 write(*,*) "7:Parallel to a bond and meantime normal to a plane defined by three atoms" 835 write(*,"(a,f8.4,a)") " 0:Set extension distance for plane type 1~5, current:",aug2D," Bohr" 836 read(*,*) plesel 837 aug2D2=aug2D !If don't draw gradient line map, needn't make the augment in the two dimension different 838 orgx2D=0D0 839 orgy2D=0D0 840 orgz2D=0D0 841 v1x=0D0 842 v1y=0D0 843 v1z=0D0 844 v2x=0D0 845 v2y=0D0 846 v2z=0D0 847 if (plesel==-10) then 848 cycle 849 else if (plesel==0) then 850 write(*,*) "Input extension distance in Bohr, e.g. 4.5" 851 read(*,*) aug2D 852 goto 401 853 else if (plesel==-1) then 854 write(*,*) "Input rotation angle of the plotting plane in degree, e.g. 30.5" 855 read(*,*) rot2Dple 856 goto 401 857 else if (plesel==1) then 858 write(*,*) "Input Z value in Bohr" 859 write(*,*) "Note: If the unit is in Angstrom, add ""a"" suffix, e.g. -1.6a" 860 read(*,*) c200tmp 861 if (index(c200tmp,'a')/=0) then 862 read(c200tmp(1:len_trim(c200tmp)-1),*) orgz2D 863 orgz2D=orgz2D/b2a 864 else 865 read(c200tmp,*) orgz2D 866 end if 867 if (ncenter==0) then 868 orgx2D=orgx 869 orgy2D=orgy 870 v1x=(endx-orgx)/(ngridnum1-1) 871 v2y=(endy-orgy)/(ngridnum2-1) 872 else 873 !Adjust aug2D/aug2D2 to make the length in two direction same 874 !Because of bug in stream() of dislin 9.5D, if the two length unequal, some region won't be shown 875 if (idrawtype==6.or.idrawtype==7) then 876 sup=(maxval(a%x)-minval(a%x))-(maxval(a%y)-minval(a%y)) 877 if (sup>0) aug2D2=aug2D+sup/2 878 if (sup<0) aug2D=aug2D-sup/2 879 end if 880 orgx2D=minval(a%x)-aug2D 881 orgy2D=minval(a%y)-aug2D2 882 v1x=(maxval(a%x)+aug2D-orgx2D)/ngridnum1 883 v2y=(maxval(a%y)+aug2D2-orgy2D)/ngridnum2 884 end if 885 else if (plesel==2) then 886 write(*,*) "Input Y value in Bohr" 887 write(*,*) "Note: If the unit is in Angstrom, add ""a"" suffix, e.g. -1.6a" 888 read(*,*) c200tmp 889 if (index(c200tmp,'a')/=0) then 890 read(c200tmp(1:len_trim(c200tmp)-1),*) orgy2D 891 orgy2D=orgy2D/b2a 892 else 893 read(c200tmp,*) orgy2D 894 end if 895 if (ncenter==0) then 896 orgx2D=orgx 897 orgz2D=orgz 898 v1x=(endx-orgx)/(ngridnum1-1) 899 v2z=(endz-orgz)/(ngridnum2-1) 900 else 901 if (idrawtype==6.or.idrawtype==7) then !adjust aug2D/aug2D2 902 sup=(maxval(a%x)-minval(a%x))-(maxval(a%z)-minval(a%z)) 903 if (sup>0) aug2D2=aug2D+sup/2 904 if (sup<0) aug2D=aug2D-sup/2 905 end if 906 orgx2D=minval(a%x)-aug2D 907 orgz2D=minval(a%z)-aug2D2 908 v1x=(maxval(a%x)+aug2D-orgx2D)/ngridnum1 909 v2z=(maxval(a%z)+aug2D2-orgz2D)/ngridnum2 910 end if 911 else if (plesel==3) then 912 write(*,*) "Input X value in Bohr" 913 write(*,*) "Note: If the unit is in Angstrom, add ""a"" suffix, e.g. -1.6a" 914 read(*,*) c200tmp 915 if (index(c200tmp,'a')/=0) then 916 read(c200tmp(1:len_trim(c200tmp)-1),*) orgx2D 917 orgx2D=orgx2D/b2a 918 else 919 read(c200tmp,*) orgx2D 920 end if 921 if (ncenter==0) then 922 orgy2D=orgy 923 orgz2D=orgz 924 v1y=(endy-orgy)/(ngridnum1-1) 925 v2z=(endz-orgz)/(ngridnum2-1) 926 else 927 if (idrawtype==6.or.idrawtype==7) then !adjust aug2D/aug2D2 928 sup=(maxval(a%y)-minval(a%y))-(maxval(a%z)-minval(a%z)) 929 if (sup>0) aug2D2=aug2D+sup/2 930 if (sup<0) aug2D=aug2D-sup/2 931 end if 932 orgy2D=minval(a%y)-aug2D 933 orgz2D=minval(a%z)-aug2D2 934 v1y=(maxval(a%y)+aug2D-orgy2D)/ngridnum1 935 v2z=(maxval(a%z)+aug2D2-orgz2D)/ngridnum2 936 end if 937 else if (plesel==4.or.plesel==5) then 938410 if (plesel==4) then 939 write(*,*) "Input the number of three atoms (e.g. 3,6,7)" 940 read(*,*) i1,i2,i3 941 if (i1==i2.or.i1==i3.or.i2==i3.or.min(i1,i2,i3)<1.or.max(i1,i2,i3)>ncenter) then 942 if (min(i1,i2,i3)<1.or.max(i1,i2,i3)>ncenter) write(*,*) "Atom indices are out of valid range, please input again" 943 if (i1==i2.or.i1==i3.or.i2==i3) write(*,*) "Atom indices are duplicated, please input again" 944 goto 410 945 end if 946 a1x=a(i1)%x 947 a1y=a(i1)%y 948 a1z=a(i1)%z 949 a2x=a(i2)%x 950 a2y=a(i2)%y 951 a2z=a(i2)%z 952 a3x=a(i3)%x 953 a3y=a(i3)%y 954 a3z=a(i3)%z 955 else if (plesel==5) then 956 write(*,*) "Input x,y,z of point 1 (in Bohr, e.g. 1.2,1.3,0.0)" 957 read(*,*) a1x,a1y,a1z 958 write(*,*) "Input x,y,z of point 2 (in Bohr)" 959 read(*,*) a2x,a2y,a2z 960 write(*,*) "Input x,y,z of point 3 (in Bohr)" 961 read(*,*) a3x,a3y,a3z 962 end if 963 v1x=a1x-a2x 964 v1y=a1y-a2y 965 v1z=a1z-a2z 966 v2x=a3x-a2x 967 v2y=a3y-a2y 968 v2z=a3z-a2z 969 rnorm1=dsqrt(v1x**2+v1y**2+v1z**2) !Norm of vector 1 970 rnorm2=dsqrt(v2x**2+v2y**2+v2z**2) 971 if (abs(v1x*v2x+v1y*v2y+v1z*v2z)/(rnorm1*rnorm2)>0.999) then 972 if (plesel==4) write(*,*) "The three atoms should not lie in the same line!" 973 if (plesel==5) write(*,*) "The three points should not lie in the same line!" 974 goto 410 975 end if 976 rangle=acos( abs(v1x*v2x+v1y*v2y+v1z*v2z)/(rnorm1*rnorm2) ) 977 if (idrawtype==6.or.idrawtype==7) then !adjust aug2D/aug2D2 978 sup=rnorm1-rnorm2*cos(pi/2-rangle) 979 if (sup>0) aug2D2=aug2D+sup/2 980 if (sup<0) aug2D=aug2D-sup/2 981 end if 982 dist1=rnorm1+2*aug2D !Total length in direction 1 983 dist2=rnorm2*cos(pi/2-rangle)+2*aug2D2 !Total length in direction 2 984! write(*,*) "angle",acos(cos(pi/2-rangle))/pi*180 985 d1=dist1/ngridnum1 !Transitional step length in direction 1 986 d2=dist2/ngridnum2 987 v1x=v1x*d1/rnorm1 !Make the norm of v1 equal to expected step lengh (d1) 988 v1y=v1y*d1/rnorm1 989 v1z=v1z*d1/rnorm1 990 schmit=(v1x*v2x+v1y*v2y+v1z*v2z)/(v1x**2+v1y**2+v1z**2) !Use schmit method to make v2 ortho to v1 991 v2x=v2x-schmit*v1x 992 v2y=v2y-schmit*v1y 993 v2z=v2z-schmit*v1z 994 rnorm2=dsqrt(v2x**2+v2y**2+v2z**2) 995 v2x=v2x*d2/rnorm2 !Make the norm of v2 equal to expected step lengh (d2) 996 v2y=v2y*d2/rnorm2 997 v2z=v2z*d2/rnorm2 998! write(*,*) "test ortho",v1x*v2x+v1y*v2y+v1z*v2z 999 orgx2D=a2x-aug2D/d1*v1x-aug2D2/d2*v2x !aug2D/d1*v1x=aug2D*(v1x/d1), v1x/d1 correspond the x component of unit vector in v1x direction 1000 orgy2D=a2y-aug2D/d1*v1y-aug2D2/d2*v2y 1001 orgz2D=a2z-aug2D/d1*v1z-aug2D2/d2*v2z 1002 else if (plesel==6) then 1003 write(*,*) "Input origin of x,y,z in Bohr, e.g. 3.5,-1,0.2" 1004 read(*,*) orgx2D,orgy2D,orgz2D 1005 write(*,*) "Input x,y,z of transitional vector 1 in Bohr, e.g. 0.08,0.03,0" 1006 read(*,*) v1x,v1y,v1z 1007 write(*,*) "Input x,y,z of transitional vector 2 in Bohr, should be orthogonal to vector 1" 1008 read(*,*) v2x,v2y,v2z 1009 d1=dsqrt(v1x**2+v1y**2+v1z**2) 1010 d2=dsqrt(v2x**2+v2y**2+v2z**2) 1011 dist1=d1*ngridnum1 1012 dist2=d2*ngridnum2 1013 a1x=orgx2D !Although a1x...a3z is no use for generate data, but these are critical for plotting atom label (subroutine drawplane) 1014 a1y=orgy2D 1015 a1z=orgz2D 1016 a2x=orgx2D+ngridnum1*v1x 1017 a2y=orgy2D+ngridnum1*v1y 1018 a2z=orgz2D+ngridnum1*v1z 1019 a3x=orgx2D+ngridnum2*v2x 1020 a3y=orgy2D+ngridnum2*v2y 1021 a3z=orgz2D+ngridnum2*v2z 1022 else if (plesel==7) then 1023 write(*,*) "Input two atoms to define the bond, e.g. 4,5" 1024 read(*,*) iatm1,iatm2 1025 write(*,*) "Input three atoms to define a plane, e.g. 1,4,7" 1026 read(*,*) jatm1,jatm2,jatm3 1027 write(*,*) "Input length of X-axis in Bohr e.g. 10" 1028 read(*,*) dist1 1029 write(*,*) "Input length of Y-axis in Bohr e.g. 8" 1030 read(*,*) dist2 1031 v1x=a(iatm2)%x-a(iatm1)%x 1032 v1y=a(iatm2)%y-a(iatm1)%y 1033 v1z=a(iatm2)%z-a(iatm1)%z 1034 call pointABCD(a(jatm1)%x,a(jatm1)%y,a(jatm1)%z,a(jatm2)%x,a(jatm2)%y,a(jatm2)%z,a(jatm3)%x,a(jatm3)%y,a(jatm3)%z,v2x,v2y,v2z,tmpD) 1035 schmit=(v1x*v2x+v1y*v2y+v1z*v2z)/(v1x**2+v1y**2+v1z**2) !Use schmit method to make v2 ortho to v1 1036 v2x=v2x-schmit*v1x 1037 v2y=v2y-schmit*v1y 1038 v2z=v2z-schmit*v1z 1039 d1=dist1/(ngridnum1-1) 1040 d2=dist2/(ngridnum2-1) 1041 rnorm1=dsqrt(v1x**2+v1y**2+v1z**2) 1042 v1x=v1x/rnorm1*d1 1043 v1y=v1y/rnorm1*d1 1044 v1z=v1z/rnorm1*d1 1045 rnorm2=dsqrt(v2x**2+v2y**2+v2z**2) 1046 v2x=v2x/rnorm2*d2 1047 v2y=v2y/rnorm2*d2 1048 v2z=v2z/rnorm2*d2 1049 orgx2D=(a(iatm1)%x+a(iatm2)%x)/2 - v1x*ngridnum1/2 - v2x*ngridnum2/2 1050 orgy2D=(a(iatm1)%y+a(iatm2)%y)/2 - v1y*ngridnum1/2 - v2y*ngridnum2/2 1051 orgz2D=(a(iatm1)%z+a(iatm2)%z)/2 - v1z*ngridnum1/2 - v2z*ngridnum2/2 1052 a1x=orgx2D !Although a1x...a3z is no use for generate data, but these are critical for plotting atom label (subroutine drawplane) 1053 a1y=orgy2D 1054 a1z=orgz2D 1055 a2x=orgx2D+ngridnum1*v1x 1056 a2y=orgy2D+ngridnum1*v1y 1057 a2z=orgz2D+ngridnum1*v1z 1058 a3x=orgx2D+ngridnum2*v2x 1059 a3y=orgy2D+ngridnum2*v2y 1060 a3z=orgz2D+ngridnum2*v2z 1061 end if 1062 1063 write(*,*) 1064 write(*,"(' X/Y/Z of origin of the plane: ',3f10.5,' Bohr')") orgx2D,orgy2D,orgz2D 1065 endx2D=orgx2D+v1x*(ngridnum1-1)+v2x*(ngridnum2-1) 1066 endy2D=orgy2D+v1y*(ngridnum1-1)+v2y*(ngridnum2-1) 1067 endz2D=orgz2D+v1z*(ngridnum1-1)+v2z*(ngridnum2-1) 1068 write(*,"(' X/Y/Z of end of the plane: ',3f10.5,' Bohr')") endx2D,endy2D,endz2D 1069 write(*,"(' X/Y/Z of translation vector 1:',3f10.5,' Bohr')") v1x,v1y,v1z 1070 write(*,"(' X/Y/Z of translation vector 2:',3f10.5,' Bohr')") v2x,v2y,v2z 1071 write(*,*) 1072 rnorm1=dsqrt(v1x**2+v1y**2+v1z**2) !The final length of vector 1 1073 rnorm2=dsqrt(v2x**2+v2y**2+v2z**2) !The final length of vector 2 1074 diff=1D-5 1075 diffv1x=diff*v1x/rnorm1 !The infinitesimal in each direction for gradient plot 1076 diffv1y=diff*v1y/rnorm1 1077 diffv1z=diff*v1z/rnorm1 1078 diffv2x=diff*v2x/rnorm2 1079 diffv2y=diff*v2y/rnorm2 1080 diffv2z=diff*v2z/rnorm2 1081 1082 !Don't directly use Multiwfn to calculate the plane data, but load from external file 1083 if (iplaneextdata==1) then 1084 open(10,file="planept.txt",status="replace") 1085 open(11,file="cubegenpt.txt",status="replace") 1086 write(10,*) ngridnum1*ngridnum2 1087 do ipt=1,ngridnum1 1088 do jpt=1,ngridnum2 1089 rnowx=orgx2D+(ipt-1)*v1x+(jpt-1)*v2x 1090 rnowy=orgy2D+(ipt-1)*v1y+(jpt-1)*v2y 1091 rnowz=orgz2D+(ipt-1)*v1z+(jpt-1)*v2z 1092 write(10,"(3f16.8)") rnowx,rnowy,rnowz 1093 write(11,"(3f16.8)") rnowx*b2a,rnowy*b2a,rnowz*b2a 1094 end do 1095 end do 1096 close(10) 1097 close(11) 1098 write(*,"(a)") " The coordinate of all points needed to be calculated have been outputted to plane.txt in current folder, the unit is in Bohr" 1099 write(*,"(a)") " cubegenpt.txt is also outputted, which is similar to plane.txt, but the unit is in Angstrom, & 1100 and there is no first line (the number of points). It can be directly utilized by cubegen" 1101 write(*,"(a)") " For example ""cubegen 0 potential CNT.fch result.cub -5 h < cubegenpt.txt""" 1102 write(*,*) 1103 write(*,"(a)") " Now input the path of the file containing function values, e.g. c:\t.txt, whose format should be identical to plane.txt, but with function value in the fourth column" 1104 write(*,"(a)") " Note: If the suffix is .cub, then the file will be recognized and loaded as output file of cubegen" 1105 do while(.true.) 1106 read(*,"(a)") c200tmp 1107 inquire(file=c200tmp,exist=alive) 1108 if (alive) exit 1109 write(*,*) "File cannot be found, input again" 1110 end do 1111 open(10,file=c200tmp,status="old") 1112 if (index(c200tmp,".cub")/=0) then !cubegen output 1113 do iskip=1,6+ncenter 1114 read(10,*) 1115 end do 1116 else 1117 read(10,*) 1118 end if 1119 do ipt=1,ngridnum1 1120 do jpt=1,ngridnum2 1121 read(10,*) rnouse,rnouse,rnouse,planemat(ipt,jpt) 1122 end do 1123 end do 1124 close(10) 1125 goto 430 1126 end if 1127 1128 !!! Start calculation of plane data now! 1129 if (infuncsel2/=4) call delvirorb(1) !Delete high-lying virtual orbitals for faster calculation, but don't do this for analyzing MO 1130 write(*,*) "Please wait..." 1131 if (infuncsel2/=12.and.expcutoff<0) write(*,"(' Note: All exponential functions exp(x) with x<',f8.3,' will be ignored ')") expcutoff 1132 CALL CPU_TIME(time_begin) 1133 call walltime(walltime1) 1134 icustom=0 1135 planemat=0D0 !For promolecular property, first clean up planemat, if not planemat may not clean 1136 if (ipromol==1) goto 421 ! To obtain promolecule property, pass the first loaded molecule 1137 !Note: If the task refers to plotting ESP map, the function of drawing gradient line will be disabled 1138 1139420 if (infuncsel2==12.and.idrawtype/=6.and.idrawtype/=7) then 1140 call planeesp(max(ngridnum1,ngridnum2)) !Special treatment to calculate ESP 1141 else if (infuncsel2==112) then !Calculate atomic Hirshfeld weighting function 1142 call genhirshplanewei(tmparrint,size(tmparrint),iHirshdenstype) 1143 ncustommap=0 1144 else if (infuncsel2==500.or.infuncsel2==510.or.infuncsel2==511.or.infuncsel2==512) then 1145 !500: Calculate rho(A)*ln[rho(A)/rho0(A)], 510: Calculate rho(A), 511: Calculate rho0(A), 512: other 1146 call genhirshplanewei(tmparrint,size(tmparrint),1) 1147 ncustommap=0 1148 do i=1,ngridnum1 !Now planemat is Hirshfeld weight of iatmentropy, and planemattmp is its density in free-state 1149 do j=1,ngridnum2 1150 rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x 1151 rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y 1152 rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z 1153 rhoA=planemat(i,j)*fdens(rnowx,rnowy,rnowz) 1154 if (infuncsel2==500) planemat(i,j)=rhoA*log(rhoA/planemattmp(i,j)) 1155 if (infuncsel2==510) planemat(i,j)=rhoA 1156 if (infuncsel2==511) planemat(i,j)=planemattmp(i,j) 1157 if (infuncsel2==512) planemat(i,j)=log(rhoA/planemattmp(i,j)) 1158 end do 1159 end do 1160 else if (infuncsel2==501) then !Calculate sum{rho(A)*ln[rho(A)/rho0(A)]} 1161 call genentroplane(1) 1162 ncustommap=0 1163 else if (infuncsel2==502) then !Calculate sum(x), where x=[rho(A)-rho0(A)]/rho(A) 1164 call genentroplane(2) 1165 ncustommap=0 1166 else if (infuncsel2==503) then !Calculate difference between total relative Shannon entropy and deformation density 1167 call genentroplane(3) 1168 ncustommap=0 1169 else 1170nthreads=getNThreads() 1171!$OMP PARALLEL DO private(i,j,rnowx,rnowy,rnowz) shared(planemat,d1add,d1min,d2add,d2min) schedule(dynamic) NUM_THREADS(nthreads) 1172 do i=1,ngridnum1 1173 do j=1,ngridnum2 1174 rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x 1175 rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y 1176 rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z 1177 if (infuncsel2==111) then 1178 planemat(i,j)=beckewei(rnowx,rnowy,rnowz,iatmbecke1,iatmbecke2) 1179 else 1180 planemat(i,j)=calcfuncall(infuncsel2,rnowx,rnowy,rnowz) 1181 if (infuncsel2==4.and.iorbsel2/=0) planemattmp(i,j)=fmo(rnowx,rnowy,rnowz,iorbsel2) !Calculate another orbital together 1182 end if 1183 if (idrawtype==6.or.idrawtype==7) then !Generate two vector to plot vector field line graph 1184 d1add(i,j)=calcfuncall(infuncsel2,rnowx+diffv1x,rnowy+diffv1y,rnowz+diffv1z) 1185 d1min(i,j)=calcfuncall(infuncsel2,rnowx-diffv1x,rnowy-diffv1y,rnowz-diffv1z) 1186 d2add(i,j)=calcfuncall(infuncsel2,rnowx+diffv2x,rnowy+diffv2y,rnowz+diffv2z) 1187 d2min(i,j)=calcfuncall(infuncsel2,rnowx-diffv2x,rnowy-diffv2y,rnowz-diffv2z) 1188 end if 1189 end do 1190 end do 1191!$OMP END PARALLEL DO 1192 end if 1193 1194421 if (ncustommap/=0) then !Calculate data for custom map 1195 if (icustom==0) then !The first time 1196 planemat_cust=planemat 1197 if (idrawtype==6.or.idrawtype==7) then 1198 d1addtmp=d1add 1199 d1mintmp=d1min 1200 d2addtmp=d2add 1201 d2mintmp=d2min 1202 end if 1203 else if (icustom/=0) then 1204 if (customop(icustom)=='+') then 1205 planemat_cust=planemat_cust+planemat 1206 if (idrawtype==6.or.idrawtype==7) then 1207 d1addtmp=d1addtmp+d1add 1208 d1mintmp=d1mintmp+d1min 1209 d2addtmp=d2addtmp+d2add 1210 d2mintmp=d2mintmp+d2min 1211 end if 1212 else if (customop(icustom)=='-') then 1213 planemat_cust=planemat_cust-planemat 1214 if (idrawtype==6.or.idrawtype==7) then 1215 d1addtmp=d1addtmp-d1add 1216 d1mintmp=d1mintmp-d1min 1217 d2addtmp=d2addtmp-d2add 1218 d2mintmp=d2mintmp-d2min 1219 end if 1220 else if (customop(icustom)=='x'.or.customop(icustom)=='*') then 1221 planemat_cust=planemat_cust*planemat 1222 if (idrawtype==6.or.idrawtype==7) then 1223 d1addtmp=d1addtmp*d1add 1224 d1mintmp=d1mintmp*d1min 1225 d2addtmp=d2addtmp*d2add 1226 d2mintmp=d2mintmp*d2min 1227 end if 1228 else if (customop(icustom)=='/') then 1229 planemat_cust=planemat_cust/planemat 1230 if (idrawtype==6.or.idrawtype==7) then 1231 d1addtmp=d1addtmp/d1add 1232 d1mintmp=d1mintmp/d1min 1233 d2addtmp=d2addtmp/d2add 1234 d2mintmp=d2mintmp/d2min 1235 end if 1236 end if 1237 end if 1238 if (icustom/=ncustommap) then !Not the final time 1239 icustom=icustom+1 1240 filename=custommapname(icustom) 1241 call dealloall 1242 write(*,"(' Loading: ',a)") trim(filename) 1243 call readinfile(filename,1) 1244 if (infuncsel2/=4) call delvirorb(0) 1245 !Generate temporary fragatm 1246 deallocate(fragatm) 1247 nfragatmnum=ncenter 1248 allocate(fragatm(nfragatmnum)) 1249 do iatm=1,ncenter 1250 fragatm(iatm)=iatm 1251 end do 1252 !Input the MO index for current file. Since the MO index may be not the same as the first loaded one 1253 if (infuncsel2==4) then 1254 write(*,"(' Input the index of the orbital to be calculated for ',a,' e.g. 3')") trim(filename) 1255 read(*,*) iorbsel 1256 end if 1257 goto 420 1258 else !The final time, reload the first loaded system 1259 planemat=planemat_cust 1260 if (idrawtype==6.or.idrawtype==7) then 1261 d1add=d1addtmp 1262 d1min=d1mintmp 1263 d2add=d2addtmp 1264 d2min=d2mintmp 1265 end if 1266 call dealloall 1267 write(*,"(' Reloading: ',a)") trim(firstfilename) 1268 call readinfile(firstfilename,1) 1269 !Recovery user-defined fragatm from the backup 1270 deallocate(fragatm) 1271 nfragatmnum=nfragatmnumbackup 1272 allocate(fragatm(nfragatmnum)) 1273 fragatm=fragatmbackup 1274 end if 1275 end if 1276 1277 if (idrawtype==6.or.idrawtype==7) then !Finish the finite differential to yield gradient 1278 gradd1=(d1add-d1min)/2/diff 1279 gradd2=(d2add-d2min)/2/diff 1280 end if 1281 CALL CPU_TIME(time_end) 1282 call walltime(walltime2) 1283 1284! open(14,file="data(h2o-O1).txt",status="old") !read data from external file, normal user don't play with this 1285! do i=1,ngridnum1 1286! do j=1,ngridnum2 1287! read(14,*) rabbish,rabbish,rabbish,planemat(i,j) 1288! end do 1289! end do 1290! close(14) 1291 1292430 if (isys==1) write(*,"(' Calculation took up CPU time',f12.2,'s, wall clock time',i8,'s')") time_end-time_begin,walltime2-walltime1 1293 write(*,*) "The minimum of data:",minval(planemat) 1294 write(*,*) "The maximum of data:",maxval(planemat) 1295 !! Set default lower and upper of Z for plot 1296 surcolorzmin=-3 1297 surcolorzmax=3 1298 if (infuncsel2==1) then 1299 drawlowlim=0.0D0 1300 drawuplim=0.65D0 1301 else if (infuncsel2==2) then 1302 drawlowlim=0.0D0 1303 drawuplim=0.65D0 1304 else if (infuncsel2==3) then 1305 drawlowlim=-8.0D0 1306 drawuplim=15.0D0 1307 else if (infuncsel2==4) then 1308 drawlowlim=-0.8D0 1309 drawuplim=0.8D0 1310 else if (infuncsel2==5) then 1311 drawlowlim=-0.1D0 1312 drawuplim=0.1D0 1313 else if (infuncsel2==8.and.ifiletype/=4) then !Nuclear ESP 1314 drawlowlim=0.0D0 1315 drawuplim=50D0 1316 surcolorzmin=0 1317 surcolorzmax=50 1318 else if (infuncsel2==8.and.ifiletype==4) then !Atomic charge ESP 1319 drawlowlim=-0.4D0 1320 drawuplim=0.4D0 1321 else if (infuncsel2==9) then 1322 drawlowlim=0.0D0 1323 drawuplim=1.0D0 1324 else if (infuncsel2==10) then 1325 drawlowlim=0.0D0 1326 drawuplim=0.8D0 1327 else if (infuncsel2==11) then 1328 drawlowlim=0.0D0 1329 drawuplim=0.1D0 1330 else if (infuncsel2==12) then 1331 drawlowlim=-0.1D0 1332 drawuplim=0.1D0 1333 else if (infuncsel2==13.or.infuncsel2==14) then 1334 drawlowlim=0D0 1335 drawuplim=1D0 1336 else if (infuncsel2==15.or.infuncsel2==16) then 1337 drawlowlim=-0.65D0 1338 drawuplim=0.65D0 1339 else if (infuncsel2==17) then 1340 drawlowlim=-0.5D0 1341 drawuplim=0.1D0 1342 else if (infuncsel2==18) then 1343 drawlowlim=0D0 1344 drawuplim=2D0 1345 else if (infuncsel2==111.or.infuncsel2==112) then !Becke/Hirshfeld weight 1346 drawlowlim=0D0 1347 drawuplim=1D0 1348 else if (infuncsel2==100.and.iuserfunc==20) then !DORI 1349 drawlowlim=0D0 1350 drawuplim=1D0 1351 else ! Include infuncsel2==100 1352 drawlowlim=0.0D0 1353 drawuplim=5.0D0 1354 end if 1355 !Set up range of X and Y axes 1356 if (plesel==1) then 1357 axlow1=orgx2D 1358 axhigh1=endx2D 1359 axlow2=orgy2D 1360 axhigh2=endy2D 1361 else if (plesel==2) then 1362 axlow1=orgx2D 1363 axhigh1=endx2D 1364 axlow2=orgz2D 1365 axhigh2=endz2D 1366 else if (plesel==3) then 1367 axlow1=orgy2D 1368 axhigh1=endy2D 1369 axlow2=orgz2D 1370 axhigh2=endz2D 1371 else if (plesel==4.or.plesel==5.or.plesel==6.or.plesel==7) then 1372 axlow1=0D0 1373 axhigh1=dist1-d1 1374 axlow2=0D0 1375 axhigh2=dist2-d2 1376 end if 1377 !Step size between labels 1378 planestpx=(axhigh1-axlow1)/7 1379 planestpy=(axhigh2-axlow2)/7 1380 planestpz=(drawuplim-drawlowlim)/10 1381 1382 XVU=150.0D0 !Reinitialize view 1383 YVU=30.0D0 1384 ZVU=7.0D0 !More suitable than 6.0D0 for drawing plane 1385 i=-1 1386 1387 idrawintbasple=0 !Refresh (3,-1) information 1388 nple3n1path=0 1389 cp2ple3n1path=0 1390 iatmcontri=0 !=0 means haven't define atomic contribution 1391 if (allocated(ple3n1path)) deallocate(ple3n1path) 1392 1393 do while(.true.) 1394 !! We have calculated planemat or x/y grad, now use the data to draw graph 1395 !! Because the max cycle is ngridnum1/2 - 1, so the upper coordinate likes dist1-d1 rather than dist1 1396 if ((i==-1.and.isilent==0).or.isavepic==1) then 1397 if (idrawtype==3.or.idrawtype==4.or.idrawtype==5) then !Draw 3D plane, first use drawplaneGUI to setup GUI, which then invokes drawplane 1398 else 1399 end if 1400 if (isavepic==1) write(*,"(a,a,a)") " Graph have been saved to ",trim(graphformat)," file with ""DISLIN"" prefix in current directory" 1401 isavepic=0 1402 end if 1403 1404 !! After show the plot once, ask user what to do next. 0 and negative options are general options 1405 write(*,*) 1406! write(*,*) "-11 Translate content of the graph" !For my personal use 1407! write(*,*) "-10 Multiply data by the data in a plane text file" 1408 if (iatmcontri==0) write(*,*) "-9 Only plot the data around certain atoms" 1409 if (iatmcontri==1) write(*,*) "-9 Recovery original plane data" 1410 if (ilenunit2D==1) write(*,*) "-8 Change length unit of the graph to Angstrom" 1411 if (ilenunit2D==2) write(*,*) "-8 Change length unit of the graph to Bohr" 1412 write(*,*) "-7 Multiply data by a factor" 1413 write(*,*) "-6 Export calculated plane data to plane.txt in current folder" 1414 write(*,*) "-5 Return to main menu" 1415 write(*,*) "-4 Switch ON/OFF of reversing ticks" 1416 write(*,"(a,i3)") " -3 Set the number of ticks between the labels, current:",iticks 1417 if (idrawtype==1.or.idrawtype==4.or.idrawtype==5) then 1418 write(*,"(a,2f7.3,f10.5)") " -2 Set stepsize in X,Y,Z axes, current:",planestpx,planestpy,planestpz 1419 else 1420 write(*,"(a,2f7.3)") " -2 Set stepsize in X and Y axes, current:",planestpx,planestpy 1421 end if 1422 write(*,*) "-1 Show the graph again" 1423 write(*,*) "0 Save the graph to a file" 1424 1425 if (idrawtype==1) then !Color-filled map 1426 if (abs(drawlowlim)<1000000.and.abs(drawuplim)<1000000) then 1427 write(*,"(a,2f15.7)") " 1 Set lower&upper limit of color scale, current:",drawlowlim,drawuplim 1428 else 1429 write(*,"(a,2(1PE15.6))") " 1 Set lower&upper limit of color scale, current:",drawlowlim,drawuplim 1430 end if 1431 if (idrawcontour==1) write(*,*) "2 Disable showing contour lines" 1432 if (idrawcontour==0) write(*,*) "2 Enable showing contour lines" 1433 write(*,*) "3 Change contour line setting" 1434 if (iatom_on_contour==0) write(*,*) "4 Enable showing atom labels and reference point" 1435 if (iatom_on_contour==1) write(*,*) "4 Disable showing atom labels and reference point" 1436 if (idrawplanevdwctr==0.and.iorbsel2==0.and.allocated(b)) write(*,*) "15 Draw a contour line of vdW surface (electron density=0.001)" !meaningless if custom operation is performed 1437 if (idrawplanevdwctr==1) write(*,*) "15 Don't draw a contour line of vdW surface" 1438 if (idrawplanevdwctr==1) write(*,*) "16 Set label size, style and color of the contour line of vdW surface" !When iorbsel2/=0, that means plot another orbital, cubmattmp will be pre-occupied 1439 if (iatom_on_contour==1) write(*,"(a,f7.3)") " 17 Set the distance criterion for showing atom labels, current:",disshowlabel 1440 if (iatmlabtype==1) write(*,*) "18 Change style of atomic labels: Only plot element symbol" 1441 if (iatmlabtype==2) write(*,*) "18 Change style of atomic labels: Only plot atomic index" 1442 if (iatmlabtype==3) write(*,*) "18 Change style of atomic labels: Plot both element symbol and atomic index" 1443 else if (idrawtype==2.or.idrawtype==6.or.idrawtype==7) then !contour map, gradient line, vector field with/without contour 1444 if (iatom_on_contour==0) write(*,*) "1 Enable showing atom labels and reference point" 1445 if (iatom_on_contour==1) write(*,*) "1 Disable showing atom labels and reference point" 1446 if (ilabel_on_contour==0) write(*,*) "2 Enable showing isovalue on contour lines" 1447 if (ilabel_on_contour==1) write(*,*) "2 Disable showing isovalue on contour lines" 1448 write(*,*) "3 Change contour line setting" 1449 if (numcp>0.or.numpath>0) write(*,*) "4 Set marks of critical points and paths" 1450 if (idrawcontour==1) write(*,*) "5 Disable showing contour lines" 1451 if (idrawcontour==0) write(*,*) "5 Enable showing contour lines" 1452 if (numcp>0.and.idrawintbasple==0) write(*,*) "6 Generate and show interbasin paths" 1453 if (numcp>0.and.idrawintbasple==1) write(*,*) "6 Delete interbasin paths" 1454 if (numcp>0.and.idrawintbasple==0) write(*,*) "7 Set stepsize and maximal iteration for interbasin path generation" 1455 1456 if (idrawtype==6) then 1457 if (igrad_arrow==0) write(*,*) "10 Show arrows" 1458 if (igrad_arrow==1) write(*,*) "10 Don't show arrows" 1459 write(*,"(a,f8.4)") " 11 Set integration step for gradient line, current:",gradplotstep 1460 write(*,"(a,f8.4)") " 12 Set interstice between gradient line, current:",gradplotdis 1461 write(*,"(a,f8.4)") " 13 Set test value for drawing a new gradient line, current:",gradplottest 1462 write(*,*) "14 Set color and line width for gradient lines" 1463 else if (idrawtype==7) then 1464 write(*,"(a,f8.4)") " 10 Set upper limit of absolute value for scaling, current:",cutgradvec 1465 if (icolorvecfield==0) write(*,*) "11 Map color to arrows" 1466 if (icolorvecfield==1) write(*,*) "11 Don't Map color to arrows" 1467 if (icolorvecfield==0) write(*,*) "12 Set color for arrow heads" !If color map was set, the color set by user themselves is nulified 1468 if (iinvgradvec==0) write(*,*) "13 Invert gradient vectors" 1469 if (iinvgradvec==1) write(*,*) "13 Don't invert gradient vectors" 1470 end if 1471 if (idrawplanevdwctr==0.and.iorbsel2==0.and.allocated(b)) write(*,*) "15 Draw a contour line of vdW surface (electron density=0.001)" !meaningless if custom operation is performed 1472 if (idrawplanevdwctr==1) write(*,*) "15 Don't draw a contour line of vdW surface" 1473 if (idrawplanevdwctr==1) write(*,*) "16 Set label size, style and color of the contour line of vdW surface" 1474 if (iatom_on_contour==1) write(*,"(a,f7.3)") " 17 Set the distance criterion for showing atom labels, current:",disshowlabel 1475 if (iatmlabtype==1) write(*,*) "18 Change style of atomic labels: Only plot element symbol" 1476 if (iatmlabtype==2) write(*,*) "18 Change style of atomic labels: Only plot atomic index" 1477 if (iatmlabtype==3) write(*,*) "18 Change style of atomic labels: Plot both element symbol and atomic index" 1478 else if (idrawtype==4.or.idrawtype==5) then !Colored relief map with/without projected color-filled map 1479 write(*,*) "1 Set color scale range for filling color" 1480 write(*,"(a,a)") " 2 Toggle drawing mesh on the surface, current: ",drawsurmesh 1481 else if (idrawtype==3) then 1482 continue 1483 end if 1484 1485 read(*,*) i 1486 1487 !Below are general options 1488 if (i==-11) then 1489 write(*,*) "Input translate extent in X and Y of the plot, respectively in Bohr" 1490 write(*,*) "e.g. 0.8,-0.2" 1491 read(*,*) transd1,transd2 1492 orgxnew=orgx2D-transd1*v1x/rnorm1-transd2*v2x/rnorm2 1493 orgynew=orgy2D-transd1*v1y/rnorm1-transd2*v2y/rnorm2 1494 orgznew=orgz2D-transd1*v1z/rnorm1-transd2*v2z/rnorm2 1495 write(*,"(a)") " In next time of plot, you should use mode 6 to define the plotting plane with below parameters (in Bohr):" 1496 write(*,"(' X/Y/Z of origin of the plane: ',3f10.5)") orgxnew,orgynew,orgznew 1497 write(*,"(' X/Y/Z of translation vector 1:',3f10.5)") v1x,v1y,v1z 1498 write(*,"(' X/Y/Z of translation vector 2:',3f10.5)") v2x,v2y,v2z 1499 else if (i==-10) then !Load plane data in another plain text file and operate to current plane data, the plane settings must be identical 1500 write(*,*) "Input file name, e.g. C:\plane.txt" 1501 read(*,"(a)") c200tmp 1502 write(*,*) "How many columns? (4 or 6. The data in the last column will be loaded)" 1503 read(*,*) ncol 1504 open(10,file=c200tmp,status="old") 1505 do i=0,ngridnum1-1 1506 do j=0,ngridnum2-1 1507 if (ncol==4) then 1508 read(10,*) tmpv,tmpv,tmpv,planemattmp(i+1,j+1) 1509 else 1510 read(10,*) tmpv,tmpv,tmpv,tmpv,tmpv,planemattmp(i+1,j+1) 1511 end if 1512 end do 1513 end do 1514 close(10) 1515 write(*,*) "Which operation? +,-,x,/" 1516 read(*,*) c200tmp(1:1) 1517 if (c200tmp(1:1)=="+") then 1518 planemat=planemat+planemattmp 1519 else if (c200tmp(1:1)=="-") then 1520 planemat=planemat-planemattmp 1521 else if (c200tmp(1:1)=="x") then 1522 planemat=planemat*planemattmp 1523 else if (c200tmp(1:1)=="/") then 1524 planemat=planemat/planemattmp 1525 end if 1526 write(*,*) "Done!" 1527 else if (i==-9) then 1528 if (iatmcontri==0) then 1529 allocate(planemat_bk(ngridnum1,ngridnum2)) 1530 planemat_bk=planemat 1531 if (allocated(tmparrint)) deallocate(tmparrint) 1532 write(*,*) "Input index of the atoms you are interested in, e.g. 2,3,7-10" 1533 read(*,"(a)") c2000tmp 1534 call str2arr(c2000tmp,ntmp) 1535 allocate(tmparrint(ntmp)) 1536 call str2arr(c2000tmp,ntmp,tmparrint) 1537 write(*,*) "Updating plane data, please wait..." 1538 do i=1,ngridnum1 !First calculate promolecular density and store it to planemat 1539 do j=1,ngridnum2 1540 rnowx=orgx2D+(i-1)*v1x+(j-1)*v2x 1541 rnowy=orgy2D+(i-1)*v1y+(j-1)*v2y 1542 rnowz=orgz2D+(i-1)*v1z+(j-1)*v2z 1543 densall=0 1544 densfrag=0 1545 do iatm=1,ncenter 1546 tmpval=calcatmdens(iatm,rnowx,rnowy,rnowz,0) 1547 densall=densall+tmpval 1548 if (any(tmparrint==iatm)) densfrag=densfrag+tmpval 1549 end do 1550 planemat(i,j)=planemat(i,j)*densfrag/densall 1551 end do 1552 end do 1553 write(*,*) "Done! The data have been updated, you can replot it" 1554 deallocate(tmparrint) 1555 iatmcontri=1 1556 else !Recovery the backed up data 1557 planemat=planemat_bk 1558 deallocate(planemat_bk) 1559 iatmcontri=0 1560 end if 1561 else if (i==-8) then 1562 if (ilenunit2D==1) then 1563 ilenunit2D=2 1564 else if (ilenunit2D==2) then 1565 ilenunit2D=1 1566 end if 1567 else if (i==-7) then 1568 write(*,*) "Input a value, e.g. 0.3" 1569 read(*,*) scaleval 1570 planemat=planemat*scaleval 1571 write(*,*) "Done!" 1572 else if (i==-6) then 1573 open(10,file="plane.txt",status="replace") 1574 do i=0,ngridnum1-1 1575 do j=0,ngridnum2-1 1576 rnowx=orgx2D+i*v1x+j*v2x 1577 rnowy=orgy2D+i*v1y+j*v2y 1578 rnowz=orgz2D+i*v1z+j*v2z 1579! if (planemat(i+1,j+1)<0.00098.or.planemat(i+1,j+1)>0.00103) cycle 1580 if (plesel==4.or.plesel==5.or.plesel==6.or.plesel==7) then 1581 write(10,"(5f10.5,1PE18.10)") rnowx*b2a,rnowy*b2a,rnowz*b2a,i*d1*b2a,j*d2*b2a,planemat(i+1,j+1) 1582 else !Plane is vertical, the coordinate in a direction is zero 1583 write(10,"(3f10.5,1PE18.10)") rnowx*b2a,rnowy*b2a,rnowz*b2a,planemat(i+1,j+1) 1584 end if 1585 end do 1586 end do 1587 close(10) 1588 if (plesel==4.or.plesel==5.or.plesel==6.or.plesel==7) then 1589 write(*,"(a)") " The column 1,2,3 correspond to actual coordinates respectively" 1590 write(*,"(a)") " The column 4,5,6 correspond to X,Y coordinates in the graph and function value respectively" 1591 else 1592 write(*,"(a)") " The column 1,2,3,4 correspond to X,Y,Z and function value respectively" 1593 end if 1594 write(*,"(a)") "Result has been ouputted to plane.txt in current folder, length unit is in Angstrom" 1595 if (idrawtype/=3.and.idrawtype/=4.and.idrawtype/=5) then 1596 if ( numcp>0.or.(numpath>0.and.imarkpath==1).or.(nple3n1path>0.and.idrawintbasple==1) ) then 1597 write(*,"(/,a)") " If also output critical points and topology/basin paths to plain text files in current folder? (y/n)" 1598 read(*,*) selectyn 1599 if (selectyn=='y'.or.selectyn=='Y') then 1600 iplaneoutall=1 !Global variable, which tells drawplane routine to export topology data 1601 iplaneoutall=0 1602 if (numcp>0) write(*,"(a)") " Critical points have been outputted to planeCP.txt in current folder. The third column is type: 1=(3,-3), 2=(3,-1), 3=(3,+1), 4=(3,+3)" 1603 if (numpath>0.and.imarkpath==1) write(*,*) "Topology paths have been outputted to planepath.txt in current folder" 1604 if (nple3n1path>0.and.idrawintbasple==1) write(*,*) "Interbasin paths have been outputted to planeinterbasin.txt in current folder" 1605 write(*,"(a)") " The first two columns correspond to X,Y coordinates in the graph, the unit is in Angstrom" 1606 end if 1607 end if 1608 end if 1609 else if (i==-5) then 1610 deallocate(planemat,planemattmp) 1611 if (allocated(planemat_bk)) deallocate(planemat_bk) 1612 idrawplanevdwctr=0 1613 iorbsel2=0 1614 exit 1615 else if (i==-4) then 1616 if (itickreverse==0) then 1617 itickreverse=1 1618 else 1619 itickreverse=0 1620 end if 1621 else if (i==-3) then 1622 write(*,*) "How many ticks between labels do you want?" 1623 read(*,*) iticks 1624 iticks=iticks+1 1625 else if (i==-2) then 1626 if (idrawtype==1.or.idrawtype==4.or.idrawtype==5) then 1627 write(*,"(a)") " Input the stepsize between the labels in X,Y,Z axes, e.g. 1.5,2.0,0.1" 1628 read(*,*) planestpx,planestpy,planestpz 1629 else 1630 write(*,"(a)") " Input the stepsize between the labels in X and Y axes, e.g. 1.5,2.0" 1631 read(*,*) planestpx,planestpy 1632 end if 1633 else if (i==-1) then 1634 cycle 1635 else if (i==0) then 1636 isavepic=1 1637 end if 1638 1639 !Shared options for idrawtype 1,2,6,7 are given first 1640 if (idrawtype==1.or.idrawtype==2.or.idrawtype==6.or.idrawtype==7) then 1641 if (i==15) then 1642 if (idrawplanevdwctr==0) then 1643 idrawplanevdwctr=1 1644 write(*,*) "Please wait..." 1645nthreads=getNThreads() 1646!$OMP PARALLEL DO private(ipt,jpt,rnowx,rnowy,rnowz) shared(planemattmp) schedule(dynamic) NUM_THREADS(nthreads) 1647 do ipt=0,ngridnum1-1 1648 do jpt=0,ngridnum2-1 1649 rnowx=orgx2D+ipt*v1x+jpt*v2x 1650 rnowy=orgy2D+ipt*v1y+jpt*v2y 1651 rnowz=orgz2D+ipt*v1z+jpt*v2z 1652 planemattmp(ipt+1,jpt+1)=fdens(rnowx,rnowy,rnowz) 1653 end do 1654 end do 1655!$OMP END PARALLEL DO 1656 write(*,*) "Done, now you can replot the graph to check effect" 1657 else if (idrawplanevdwctr==1) then 1658 idrawplanevdwctr=0 1659 end if 1660 else if (i==16) then 1661 write(*,*) "Input the size of the label, e.g. 30" 1662 write(*,*) "Note: If you input 0, then the label will not be shown" 1663 read(*,*) ivdwctrlabsize 1664 write(*,*) "Select color of the contour line:" 1665 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 1666 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 1667 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 1668 read(*,*) ivdwclrindctr 1669 write(*,*) "Input the width of the contour line, e.g. 10" 1670 read(*,*) iwidthvdwctr 1671 write(*,*) "Input length of line segment and interstice" 1672 write(*,*) "e.g. 1,0 means solid line; 1,10 means DOT; 10,15 means DASH" 1673 write(*,*) " 10,25 means DASH with larger interstice" 1674 read(*,*) vdwctrstyle 1675 write(*,*) "Done, now you can replot the graph to check effect" 1676 else if (i==17) then 1677 write(*,"(a)") " Note: If the distance between an atom nucleus/critical point and the plane of interest is less than this criterion, & 1678 the atom/critical point will be labelled on the graph" 1679 write(*,*) "Input the criterion value in Bohr, e.g. 0.5" 1680 read(*,*) disshowlabel 1681 write(*,"(a)") " If also show the labels of the atoms that beyond this criterion as light face type? (y/n)" 1682 read(*,*) selectyn 1683 if (selectyn=='y'.or.selectyn=='Y') then 1684 iatom_on_contour_far=1 1685 else 1686 iatom_on_contour_far=0 1687 end if 1688 else if (i==18) then 1689 write(*,*) "1: Only plot element symbol" 1690 write(*,*) "2: Only plot atomic index" 1691 write(*,*) "3: Plot both element symbol and atomic index" 1692 write(*,*) "Note that the default value can be set by ""iatmlabtype"" in settings.ini" 1693 read(*,*) iatmlabtype 1694 end if 1695 1696 !Options only for idrawtype 1===================== 1697 if (idrawtype==1) then 1698 if (i==1) then 1699 write(*,*) "Input lower & upper limit of Z (e.g. -0.3,0.3)" 1700 read(*,*) drawlowlim,drawuplim 1701 else if (i==2) then 1702 if (idrawcontour==1) then 1703 idrawcontour=0 1704 else if (idrawcontour==0) then 1705 idrawcontour=1 1706 end if 1707 else if (i==3) then !! Change isovalues of contour line 1708 call setctr 1709 else if (i==4) then 1710 if (iatom_on_contour==1) then 1711 iatom_on_contour=0 1712 else if (iatom_on_contour==0) then 1713 iatom_on_contour=1 1714 write(*,*) "Use which color for labelling atoms?" 1715 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 1716 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 1717 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 1718 read(*,*) iclrindatmlab 1719 end if 1720 end if 1721 !Options only for idrawtype 2,6,7======================== 1722 else if (idrawtype==2.or.idrawtype==6.or.idrawtype==7) then 1723 !General option for idrawtype 2,6,7 1724 if (i==1) then 1725 if (iatom_on_contour==1) then 1726 iatom_on_contour=0 1727 else if (iatom_on_contour==0) then 1728 iatom_on_contour=1 1729 write(*,*) "Use which color?" 1730 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 1731 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 1732 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 1733 read(*,*) iclrindatmlab 1734 end if 1735 else if (i==2) then 1736 if (ilabel_on_contour==1) then 1737 ilabel_on_contour=0 1738 else if (ilabel_on_contour==0) then 1739 ilabel_on_contour=1 1740 write(*,*) "Input label size e.g. 30" 1741 read(*,*) ictrlabsize 1742 write(*,"(a)") " Hint: The number of digits after the decimal point of label on contour lines can be set by ""numdigctr"" in settings.ini" 1743 end if 1744 else if (i==3) then !! Change isovalues of contour line 1745 call setctr 1746 else if (i==4) then 1747 call settopomark 1748 else if (i==5) then 1749 if (idrawcontour==1) then 1750 idrawcontour=0 1751 else if (idrawcontour==0) then 1752 idrawcontour=1 1753 end if 1754 else if (i==6) then !Generate paths from (3,-1) 1755 if (idrawintbasple==0) then 1756 do icp=1,numcp 1757 if (CPtype(icp)==2) then 1758 cpx=CPpos(1,icp) 1759 cpy=CPpos(2,icp) 1760 cpz=CPpos(3,icp) 1761 if (plesel==1) then 1762 if (abs(cpz-orgz2D) > disshowlabel) cycle 1763 else if (plesel==2) then 1764 if (abs(cpy-orgy2D) > disshowlabel) cycle 1765 else if (plesel==3) then 1766 if (abs(cpx-orgx2D) > disshowlabel) cycle 1767 else if (plesel==4.or.plesel==5.or.plesel==6.or.plesel==7) then 1768 call pointprjple(a1x,a1y,a1z,a2x,a2y,a2z,a3x,a3y,a3z,cpx,cpy,cpz,prjx,prjy,prjz) 1769 if ( (cpx-prjx)**2+(cpy-prjy)**2+(cpz-prjz)**2 > disshowlabel**2) cycle 1770 end if 1771 nple3n1path=nple3n1path+1 1772 cp2ple3n1path(icp)=nple3n1path !Default is zero, means this CP is not on the given plane and has no corresponding interbasin path 1773 end if 1774 end do 1775 if (nple3n1path>0) then 1776 idrawintbasple=1 1777 write(*,"(' Found',i8,' (3,-1) CPs in the plane')") nple3n1path 1778 allocate(ple3n1path(3,n3n1plept,2,nple3n1path)) 1779 write(*,*) "Generating interbasin paths from (3,-1) CPs, Please wait..." 1780nthreads=getNThreads() 1781!$OMP PARALLEL DO SHARED(numcp) PRIVATE(icp) schedule(dynamic) NUM_THREADS(nthreads) 1782 do icp=1,numcp 1783 if (cp2ple3n1path(icp)/=0) call gen3n1plepath(ifunctopo,icp,cp2ple3n1path(icp)) 1784 ! write(*,"('Finished the in-plane path generation from (3,-1)',i8)") icp 1785 end do 1786!$OMP END PARALLEL DO 1787 else 1788 write(*,*) "No (3,-1) CP is closed to the plane" 1789 end if 1790 else if (idrawintbasple==1) then 1791 idrawintbasple=0 1792 nple3n1path=0 1793 cp2ple3n1path=0 1794 if (allocated(ple3n1path)) deallocate(ple3n1path) 1795 end if 1796 else if (i==7) then 1797 write(*,*) "Input stepsize (Bohr) and maximal iterations" 1798 write(*,"(a,f8.5,',',i6)") " Note: Current values:",ple3n1pathstpsiz,n3n1plept 1799 read(*,*) ple3n1pathstpsiz,n3n1plept 1800 end if 1801 1802 !Option only for idrawtype 6 1803 if (idrawtype==6) then 1804 if (i==10) then 1805 if (igrad_arrow==1) then 1806 igrad_arrow=0 1807 else if (igrad_arrow==0) then 1808 igrad_arrow=1 1809 end if 1810 else if (i==11) then 1811 write(*,*) "Input a value, default value is 0.002D0" 1812 read(*,*) gradplotstep 1813 else if (i==12) then 1814 write(*,"(a,f8.4)") "Input a value, should be larger than",gradplotstep 1815 read(*,*) gradplotdis 1816 else if (i==13) then 1817 write(*,*) "Input a value, default value is 0.2" 1818 read(*,*) gradplottest 1819 else if (i==14) then 1820 write(*,*) "Use which color?" 1821 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 1822 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 1823 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 1824 read(*,*) iclrindgradline 1825 write(*,*) "Input line width, e.g. 5 (default value is 1)" 1826 read(*,*) iwidthgradline 1827 end if 1828 !Individual option for idrawtype 7 1829 else if (idrawtype==7) then 1830 if (i==10) then 1831 write(*,*) "Input a value" 1832 read(*,*) cutgradvec 1833 else if (i==11) then 1834 if (icolorvecfield==1) then 1835 icolorvecfield=0 1836 else if (icolorvecfield==0) then 1837 icolorvecfield=1 1838 end if 1839 else if (i==12) then 1840 write(*,*) "Input color index, e.g. 1 =black, 50 =blue, 150 =green, 250 =red" 1841 read(*,*) vecclrind 1842 else if (i==13) then 1843 if (iinvgradvec==1) then 1844 iinvgradvec=0 1845 else if (iinvgradvec==0) then 1846 iinvgradvec=1 1847 end if 1848 end if 1849 end if 1850 end if 1851 1852 !Options for idrawtype 4,5================= 1853 else if (idrawtype==4.or.idrawtype==5) then 1854 if (i==1) then 1855 write(*,*) "Input lower & upper limits of color scale for shading surface, e.g. -0.1,0.3" 1856 write(*,"(a,f14.6,a,f14.6)") " Present value: from",surcolorzmin," to",surcolorzmax 1857 read(*,*) surcolorzmin,surcolorzmax 1858 if (idrawtype==5) then 1859 write(*,*) "Input lower & upper limits of color scale for projected map, e.g. -0.1,0.3" 1860 write(*,"(a,f14.6,a,f14.6)") " Present value: from",drawlowlim," to",drawuplim 1861 read(*,*) drawlowlim,drawuplim 1862 end if 1863 else if (i==2) then 1864 if (drawsurmesh=="ON ") then 1865 drawsurmesh="OFF" 1866 else 1867 drawsurmesh="ON" 1868 end if 1869 end if 1870 else if (idrawtype==3) then !No options for map type 3 currently 1871 continue 1872 end if 1873 end do 1874 1875 1876 1877!!!-------------------------------------------------------- 1878!!!-------------------------------------------------------- 1879!5!------------------- Calculate, show and output grid file 1880!!!-------------------------------------------------------- 1881!!!-------------------------------------------------------- 1882 1883else if (infuncsel1==5) then 1884 ncustommap=0 !Clean custom operation setting that possibly defined by other modules 1885 if (allocated(custommapname)) deallocate(custommapname) 1886 if (allocated(customop)) deallocate(customop) 1887 write(*,*) "-10 Return to main menu" 1888 write(*,*) "-2 Obtain of deformation property" 1889 write(*,*) "-1 Obtain of promolecule property" 1890 write(*,*) "0 Set custom operation" 1891500 call selfunc_interface(infuncsel2) 1892 1893 if (infuncsel2==0.or.infuncsel2==-1.or.infuncsel2==-2) then 1894 if (infuncsel2==0) call customplotsetup 1895 if (infuncsel2==-1) then 1896 ipromol=1 1897 else 1898 ipromol=0 1899 end if 1900 if (infuncsel2==-1.or.infuncsel2==-2) call setPromol 1901 write(*,*) "-10 Return to main menu" 1902 goto 500 1903 else if (infuncsel2==-10) then 1904 cycle 1905 else if (infuncsel2==111) then !Calculate Becke weighting function 1906 write(*,*) "Input indices of two atoms to calculate Becke overlap weight, e.g. 1,4" 1907 write(*,*) "or input index of an atom and zero to calculate Becke atomic weight, e.g. 5,0" 1908 read(*,*) iatmbecke1,iatmbecke2 1909 else if (infuncsel2==112) then !Calculate Hirshfeld weighting function 1910 if (allocated(tmparrint)) deallocate(tmparrint) 1911 write(*,*) "Input index of the atoms you are interested in, e.g. 2,3,7-10" 1912 read(*,"(a)") c2000tmp 1913 call str2arr(c2000tmp,ntmp) 1914 allocate(tmparrint(ntmp)) 1915 call str2arr(c2000tmp,ntmp,tmparrint) 1916 write(*,"(a)") " How to generate the atomic densities that used in the calculation of Hirshfeld weight?" 1917 write(*,*) "1 Based on atomic .wfn files" 1918 write(*,*) "2 Based on built-in atomic densities (see Appendix 3 of the manual for detail)" 1919 read(*,*) iHirshdenstype 1920 if (iHirshdenstype==1) call setpromol 1921 end if 1922 1923 call setgrid(1,igridsel) 1924 1925 if (igridsel==100) then !Calculate value on a set of points loaded from external file 1926 if (allocated(extpttmp)) deallocate(extpttmp) 1927 if (ncustommap/=0) allocate(extpttmp(numextpt)) !!! temp file for difference cube 1928 if (infuncsel2/=4) call delvirorb(1) !Delete high-lying virtual orbitals for faster calculation 1929 icustom=0 1930 extpt(:,4)=0D0 1931 call walltime(walltime1) 1932 CALL CPU_TIME(time_begin) 1933 if (ipromol==1) goto 509 !Calculate promolecular property, so skip the first time calculation (namely for the whole system) 1934 508 continue 1935nthreads=getNThreads() 1936!$OMP PARALLEL DO SHARED(extpt) PRIVATE(iextpt) schedule(dynamic) NUM_THREADS(nthreads) 1937 do iextpt=1,numextpt !Calculate function value 1938 extpt(iextpt,4)=calcfuncall(infuncsel2,extpt(iextpt,1),extpt(iextpt,2),extpt(iextpt,3)) 1939 end do 1940!$OMP END PARALLEL DO 1941 509 if (ncustommap/=0) then !cycling will stop when all the file have been dealed 1942 if (icustom==0) then 1943 !Note: For promolecular property, x,y,z hasn't been saved in extpt at first time, while after calculation of atoms, extpt already has %x,%y,%z 1944 extpttmp(:)=extpt(:,4) !first time 1945 else if (icustom/=0) then !not first time 1946 if (customop(icustom)=='+') extpttmp(:)=extpttmp(:)+extpt(:,4) 1947 if (customop(icustom)=='-') extpttmp(:)=extpttmp(:)-extpt(:,4) 1948 if (customop(icustom)=='x'.or.customop(icustom)=='*') extpttmp(:)=extpttmp(:)*extpt(:,4) 1949 if (customop(icustom)=='/') extpttmp(:)=extpttmp(:)/extpt(:,4) 1950 end if 1951 if (icustom/=ncustommap) then 1952 icustom=icustom+1 1953 filename=custommapname(icustom) 1954 call dealloall 1955 write(*,"(' Loading: ',a)") trim(filename) 1956 call readinfile(filename,1) 1957 if (infuncsel2/=4) call delvirorb(0) 1958 !Generate temporary fragatm 1959 deallocate(fragatm) 1960 nfragatmnum=ncenter 1961 allocate(fragatm(nfragatmnum)) 1962 do iatm=1,ncenter 1963 fragatm(iatm)=iatm 1964 end do 1965 !Input the MO index for current file. Since the MO index may be not the same as the first loaded one 1966 if (infuncsel2==4) then 1967 write(*,"(' Input the index of the orbital to be calculated for ',a,' e.g. 3')") trim(filename) 1968 read(*,*) iorbsel 1969 end if 1970 goto 508 1971 else if (icustom==ncustommap) then !last time 1972 extpt(:,4)=extpttmp(:) 1973 call dealloall 1974 write(*,"(' Reloading: ',a)") trim(firstfilename) 1975 call readinfile(firstfilename,1) 1976 !Recovery user-defined fragatm from the backup 1977 deallocate(fragatm) 1978 nfragatmnum=nfragatmnumbackup 1979 allocate(fragatm(nfragatmnum)) 1980 fragatm=fragatmbackup 1981 end if 1982 end if 1983 CALL CPU_TIME(time_end) 1984 call walltime(walltime2) 1985 if (isys==1) write(*,"(' Calculation is finished, took up CPU time',f10.2,'s, wall clock time',i9,'s')") time_end-time_begin,walltime2-walltime1 1986 write(*,"(a)") " Output the points with function values to which file? e.g. c:\ltwd.txt" 1987 read(*,"(a)") c200tmp 1988 open(10,file=c200tmp,status="replace") 1989 write(10,"(i10)") numextpt 1990 do iextpt=1,numextpt 1991 write(10,"(3f13.7,E20.10)") extpt(iextpt,1:3),extpt(iextpt,4) 1992 end do 1993 close(10) 1994 write(*,"(a)") " Done! In this file the first line is the number of points, Column 1~4 correspond to X,Y,Z coordinates and function values, respectively. All units are in a.u." 1995 1996 else !Calculate grid data 1997 if (allocated(cubmat)) deallocate(cubmat) 1998 allocate(cubmat(nx,ny,nz)) 1999 if (allocated(cubmattmp)) deallocate(cubmattmp) 2000 if (ncustommap/=0) allocate(cubmattmp(nx,ny,nz)) !!! temp file for difference cube 2001 if (infuncsel2/=4) call delvirorb(1) !Delete high-lying virtual orbitals for faster calculation 2002 !!!!! Calculate grid data 2003 if (infuncsel2==111) then !Becke's weight 2004 do k=1,nz 2005 tmpz=orgz+(k-1)*dz 2006 do j=1,ny 2007 tmpy=orgy+(j-1)*dy 2008 do i=1,nx 2009 tmpx=orgx+(i-1)*dx 2010 cubmat(i,j,k)=beckewei(tmpx,tmpy,tmpz,iatmbecke1,iatmbecke2) 2011 end do 2012 end do 2013 end do 2014 else if (infuncsel2==112) then !Hirshfeld weight 2015 ncustommap=0 2016 call genhirshcubewei(tmparrint,size(tmparrint),iHirshdenstype) 2017 else if (infuncsel2==120) then !Calculate and output three components of Steric force to plain text file 2018 open(20,file="stericforce.txt",status="replace") 2019 do k=1,nz 2020 write(*,"(' Finished:',i5,' /',i5)") k,nz 2021 tmpz=orgz+(k-1)*dz 2022 do j=1,ny 2023 tmpy=orgy+(j-1)*dy 2024 do i=1,nx 2025 tmpx=orgx+(i-1)*dx 2026 call stericderv(tmpx,tmpy,tmpz,tmpvec) 2027 write(20,"(7f12.6)") (orgx+(i-1)*dx)*b2a,(orgy+(j-1)*dy)*b2a,(orgz+(k-1)*dz)*b2a,-tmpvec, dsqrt(sum(tmpvec**2)) 2028 end do 2029 end do 2030 end do 2031 close(20) 2032 write(*,*) "Done, the results have been outputted to stericforce.txt in current folder" 2033 write(*,"(a)") " Columns 1,2,3 correspond to X,Y,Z coordinates, 4,5,6 correspond to steric force component in X,Y,Z. The last column denotes magnitude of steric force" 2034 write(*,*) 2035 read(*,*) 2036 else !Common cases 2037 cubmat=0D0 2038 icustom=0 2039 if (ipromol==1) goto 511 !Calculate promolecular property, so skip the first time calculation (namely for the whole system) 2040 510 call savecubmat(infuncsel2,0,iorbsel) !Save data to cubmat matrix 2041 511 if (ncustommap/=0) then !Calculate data for custom cube, cycling stop when all the file have been dealed 2042 if (icustom==0) then 2043 !Note: For promolecular property, x,y,z hasn't been saved in cubmat at first time, while after calculation of atoms, cubmat already has %x,%y,%z 2044 cubmattmp=cubmat !first time 2045 else if (icustom/=0) then !not first time 2046 if (customop(icustom)=='+') cubmattmp=cubmattmp+cubmat 2047 if (customop(icustom)=='-') cubmattmp=cubmattmp-cubmat 2048 if (customop(icustom)=='x'.or.customop(icustom)=='*') cubmattmp=cubmattmp*cubmat 2049 if (customop(icustom)=='/') cubmattmp=cubmattmp/cubmat 2050 end if 2051 if (icustom/=ncustommap) then 2052 icustom=icustom+1 2053 filename=custommapname(icustom) 2054 call dealloall 2055 write(*,"(' Loading: ',a)") trim(filename) 2056 call readinfile(filename,1) 2057 if (infuncsel2/=4) call delvirorb(0) 2058 !Generate temporary fragatm 2059 deallocate(fragatm) 2060 nfragatmnum=ncenter 2061 allocate(fragatm(nfragatmnum)) 2062 do iatm=1,ncenter 2063 fragatm(iatm)=iatm 2064 end do 2065 !Input the MO index for current file. Since the MO index may be not the same as the first loaded one 2066 if (infuncsel2==4) then 2067 write(*,"(' Input the index of the orbital to be calculated for ',a,' e.g. 3')") trim(filename) 2068 read(*,*) iorbsel 2069 end if 2070 goto 510 2071 else if (icustom==ncustommap) then !last time 2072 cubmat=cubmattmp 2073 call dealloall 2074 write(*,"(' Reloading: ',a)") trim(firstfilename) 2075 call readinfile(firstfilename,1) 2076 !Recovery user-defined fragatm from the backup 2077 deallocate(fragatm) 2078 nfragatmnum=nfragatmnumbackup 2079 allocate(fragatm(nfragatmnum)) 2080 fragatm=fragatmbackup 2081 end if 2082 end if 2083 end if 2084 2085 !! Output result 2086 outcubfile="griddata.cub" !General name 2087 if (infuncsel2==1) then 2088 outcubfile="density.cub" 2089 dipx=0 2090 dipy=0 2091 dipz=0 2092 do k=1,nz 2093 do j=1,ny 2094 do i=1,nx 2095 dipx=dipx-cubmat(i,j,k)*(orgx+(i-1)*dx) 2096 dipy=dipy-cubmat(i,j,k)*(orgy+(j-1)*dy) 2097 dipz=dipz-cubmat(i,j,k)*(orgz+(k-1)*dz) 2098 end do 2099 end do 2100 end do 2101 dipx=sum(a%charge*a%x)+dipx*dx*dy*dz 2102 dipy=sum(a%charge*a%y)+dipy*dx*dy*dz 2103 dipz=sum(a%charge*a%z)+dipz*dx*dy*dz 2104 write(*,*) 2105 write(*,*) "System dipole moment in a.u. (e/Bohr) and Debye, respectively:" 2106 write(*,"(' X component is',2f12.6)") dipx,dipx*8.47835281D-30*2.99792458D+29 2107 write(*,"(' Y component is',2f12.6)") dipy,dipy*8.47835281D-30*2.99792458D+29 2108 write(*,"(' Z component is',2f12.6)") dipz,dipz*8.47835281D-30*2.99792458D+29 2109 write(*,"(' Total magnitude is ',2f12.6)") dsqrt(dipx**2+dipy**2+dipz**2),dsqrt(dipx**2+dipy**2+dipz**2)*8.47835281D-30*2.99792458D+29 2110 write(*,*) 2111 else if (infuncsel2==2) then 2112 outcubfile="gradient.cub" 2113 else if (infuncsel2==3) then 2114 outcubfile="laplacian.cub" 2115 else if (infuncsel2==4) then 2116 outcubfile="MOvalue.cub" 2117 else if (infuncsel2==5) then 2118 outcubfile="spindensity.cub" 2119 else if (infuncsel2==6) then 2120 outcubfile="K(r).cub" 2121 else if (infuncsel2==7) then 2122 outcubfile="G(r).cub" 2123 else if (infuncsel2==8) then 2124 outcubfile="nucleiesp.cub" 2125 else if (infuncsel2==9) then 2126 outcubfile="ELF.cub" 2127 sur_value=0.7 2128 else if (infuncsel2==10) then 2129 outcubfile="LOL.cub" 2130 sur_value=0.5 2131 else if (infuncsel2==11) then 2132 outcubfile="infoentro.cub" 2133 else if (infuncsel2==12) then 2134 outcubfile="totesp.cub" 2135 else if (infuncsel2==13) then 2136 sur_value=0.5 2137 outcubfile="RDG.cub" 2138 else if (infuncsel2==14) then 2139 sur_value=0.4 2140 outcubfile="RDGprodens.cub" 2141 else if (infuncsel2==15) then 2142 outcubfile="signlambda2rho.cub" 2143 else if (infuncsel2==16) then 2144 outcubfile="signlambda2rhoprodens.cub" 2145 else if (infuncsel2==17) then 2146 outcubfile="fermihole.cub" 2147 else if (infuncsel2==18) then 2148 outcubfile="avglocion.cub" 2149 else if (infuncsel2==19) then 2150 outcubfile="srcfunc.cub" 2151 else if (infuncsel2==20) then 2152 outcubfile="EDR.cub" 2153 else if (infuncsel2==21) then 2154 outcubfile="EDRDmax.cub" 2155 else if (infuncsel2==100) then 2156 outcubfile="userfunc.cub" 2157 else if (infuncsel2==111) then 2158 outcubfile="Becke.cub" 2159 else if (infuncsel2==112) then 2160 outcubfile="Hirshfeld.cub" 2161 end if 2162 2163 temp=minval(cubmat) 2164 call findvalincub(cubmat,temp,i,j,k) 2165 write(*,"(' The minimum is',D16.8,' at',3f10.5,' (Bohr)')") temp,orgx+(i-1)*dx,orgy+(j-1)*dy,orgz+(k-1)*dz 2166 temp=maxval(cubmat) 2167 call findvalincub(cubmat,temp,i,j,k) 2168 write(*,"(' The maximum is',D16.8,' at',3f10.5,' (Bohr)')") temp,orgx+(i-1)*dx,orgy+(j-1)*dy,orgz+(k-1)*dz 2169 write(*,"(' Summing up all value and multiply differential element:')") 2170 write(*,*) sum(cubmat)*dx*dy*dz 2171 write(*,"(' Summing up positive value and multiply differential element:')") 2172 write(*,*) sum(cubmat,mask=cubmat>0)*dx*dy*dz 2173 write(*,"(' Summing up negative value and multiply differential element:')") 2174 write(*,*) sum(cubmat,mask=cubmat<0)*dx*dy*dz 2175 2176 !Reinitialize plot parameter 2177 bondcrit=1.15D0 2178 textheigh=38.0D0 2179 ratioatmsphere=1.0D0 2180 bondradius=0.2D0 2181 ishowatmlab=1 2182 ishowaxis=1 2183 idrawmol=1 2184 isosurshowboth=1 2185 ishowdatarange=0 2186 2187 do while(.true.) 2188 write(*,*) 2189 write(*,*) "-1 Show isosurface graph" 2190 write(*,*) "0 Return to main menu" 2191 write(*,*) "1 Save graph of isosurface to file in current folder" 2192 write(*,*) "2 Export data to Gaussian cube file in current folder" 2193 write(*,*) "3 Export data to formatted text file in current folder" 2194 write(*,"(a,f10.5)") " 4 Set the value of isosurface to be shown, current:",sur_value 2195 write(*,*) "5 Multiply all grid data by a factor" 2196 write(*,*) "6 Divide all grid data by a factor" 2197 write(*,*) "7 Add a value to all grid data" 2198 write(*,*) "8 Substract a value from all grid data" 2199 read(*,*) i 2200 2201 if (i==-1) then 2202 else if (i==0) then 2203 exit 2204 else if (i==1) then 2205 idrawisosur=1 2206 isavepic=1 2207 isavepic=0 2208 write(*,*) "Graph has been saved to current folder with ""DISLIN"" prefix" 2209 else if (i==2) then 2210 open(10,file=outcubfile,status="replace") 2211 call outcube(cubmat,nx,ny,nz,orgx,orgy,orgz,dx,dy,dz,10) 2212 close(10) 2213 write(*,"(' New grid file has been outputted to: ',a35)") adjustl(outcubfile) 2214 else if (i==3) then 2215 open(10,file="output.txt",status="replace") 2216 write(*,*) "Outputting output.txt in current directory..." 2217 do i=1,nx 2218 do j=1,ny 2219 do k=1,nz 2220 write(10,"(3f12.6,2x,1PD15.8)") (orgx+(i-1)*dx)*b2a,(orgy+(j-1)*dy)*b2a,(orgz+(k-1)*dz)*b2a,cubmat(i,j,k) 2221 end do 2222 end do 2223 end do 2224 close(10) 2225 write(*,*) "Output finished, column 1/2/3/4 correspond to X/Y/Z/Value, unit is Angstrom" 2226 else if (i==4) then 2227 write(*,*) "Input the value of isosurface" 2228 read(*,*) sur_value 2229 else if (i==5) then 2230 write(*,*) "Input a value" 2231 read(*,*) tmpval 2232 cubmat=cubmat*tmpval 2233 else if (i==6) then 2234 write(*,*) "Input a value" 2235 read(*,*) tmpval 2236 cubmat=cubmat/tmpval 2237 else if (i==7) then 2238 write(*,*) "Input a value" 2239 read(*,*) tmpval 2240 cubmat=cubmat+tmpval 2241 else if (i==8) then 2242 write(*,*) "Input a value" 2243 read(*,*) tmpval 2244 cubmat=cubmat-tmpval 2245 end if 2246 end do 2247 end if 2248 2249!!!--------------------------------------- 2250!6!!------------------- Check & Modify wavefunction or show GTF/Orbital information 2251else if (infuncsel1==6) then 2252 call modwfn 2253 2254 2255!!!--------------------------------------- 2256!7!!------------------- Population analysis 2257else if (infuncsel1==7) then 2258 call popana 2259 2260 2261!!!--------------------------------------- 2262!8!!------------------- Orbital composition analysis 2263else if (infuncsel1==8) then 2264 call compana 2265 2266 2267!!!--------------------------------------- 2268!9!!------------------- Bond order analysis 2269else if (infuncsel1==9) then 2270 call bondana 2271 2272 2273!!!--------------------------------------- 2274!10!!------------------- Plot DOS 2275else if (infuncsel1==10) then 2276 if (ifiletype/=0.or.ifiletype/=1.or.ifiletype/=9) then 2277 call plotdos 2278 else 2279 write(*,"(a,/)") " Error: This function is only available for input file containing basis function information & 2280 (e.g. .fch/.molden/.gms) and plain text file with energy levels!" 2281 end if 2282 2283 2284!!!--------------------------------------- 2285!11!!------------------- Plot spectrums 2286else if (infuncsel1==11) then 2287 call plotspectrum 2288 2289 2290!!!--------------------------------------- 2291!12!!------------------- Molecular surface analysis 2292else if (infuncsel1==12) then 2293 call surfana 2294 2295 2296!!!--------------------------------------- 2297!13!!------------------- Process grid data 2298else if (infuncsel1==13) then 2299 if (.not.allocated(cubmat)) then 2300 write(*,"(a)") " Error: Grid data was not loaded or generated! If you want to load a grid data now, input its path, e.g. C:\nico.cub, else input 0 to exit" 2301 do while(.true.) 2302 read(*,"(a)") c200tmp 2303 if (c200tmp(1:1)=='0') then 2304 exit 2305 else 2306 inquire(file=c200tmp,exist=alive) 2307 if (alive) then 2308 inamelen=len_trim(c200tmp) 2309 !Only load grid data, do not pertube other variables 2310 if (c200tmp(inamelen-2:inamelen)=="grd") then 2311 call readgrd(c200tmp,1,1) 2312 else if (c200tmp(inamelen-2:inamelen)=="cub".or.c200tmp(inamelen-3:inamelen)=="cube") then 2313 call readcube(c200tmp,1,1) 2314 else 2315 write(*,*) "Error: Unknown file type, input again" 2316 cycle 2317 end if 2318 call procgriddata 2319 exit 2320 else 2321 write(*,*) "Cannot find the file, input again" 2322 end if 2323 end if 2324 end do 2325 else 2326 call procgriddata 2327 end if 2328 2329 2330!!!--------------------------------------- 2331!14!!------------------- Adaptive natural density partitioning (AdNDP) 2332else if (infuncsel1==14) then 2333 call AdNDP 2334 2335 2336!!!--------------------------------------- 2337!15!!------------------- Integrate fuzzy atomic space 2338else if (infuncsel1==15) then 2339 call intatomspace(0) 2340else if (infuncsel1==-15) then 2341 call fuzzySBL 2342 2343!!!--------------------------------------- 2344!16!!------------------- Charge decomposition analysis 2345else if (infuncsel1==16) then 2346 write(*,*) "Citation of original GCDA and CDA used in Multiwfn, respectively:" 2347 write(*,"(a)") " GCDA: Meng Xiao, Tian Lu, Generalized Charge Decomposition Analysis (GCDA) Method, J. Adv. Phys. Chem., 4, 111-124 (2015), http://dx.doi.org/10.12677/JAPC.2015.44013" 2348 write(*,"(a)") " CDA: Stefan Dapprich, Gernot Frenking, J. Phys. Chem., 99, 9352-9362 (1995)" 2349 write(*,*) 2350 call CDA 2351 2352 2353!!!--------------------------------------- 2354!17!!------------------- Basin integration 2355else if (infuncsel1==17) then 2356 call basinana 2357 2358 2359!!!--------------------------------------- 2360!18!!------------------- Electron excitation analysis 2361else if (infuncsel1==18) then 2362 do while(.true.) 2363 write(*,*) " ------------ Electron excitation analyses ------------ " 2364 write(*,*) "0 Return" 2365 write(*,"(a)") " 1 Analyze and visualize hole-electron distribution, transition dipole moment and transition density" 2366 write(*,*) "2 Plot transition density matrix as color-filled map" 2367 write(*,*) "3 Analyze charge-transfer based on density difference grid data (JCTC,7,2498)" 2368 write(*,*) "4 Calculate delta_r index to measure charge-transfer length (JCTC,9,3118)" 2369 write(*,*) "5 Calculate transition dipole moments between all excited states" 2370 write(*,*) "6 Generate natural transition orbitals (NTOs)" 2371 write(*,*) "7 Calculate ghost-hunter index (JCC,38,2151)" 2372 write(*,*) "8 Calculate interfragment charge transfer in electronic excitation" 2373 2374 read(*,*) isel 2375 if (isel==0) then 2376 goto 10 2377 else if (isel==1) then 2378 call hetransdipdens(1) 2379 else if (isel==2) then 2380 call plottransdensmat 2381 else if (isel==3) then 2382 if (allocated(cubmat)) then 2383 call CTanalyze 2384 else 2385 write(*,"(a,/)") " Error: Grid data of electron density difference must be calculated by main function 5 or loaded from external file first!" 2386 end if 2387 else if (isel==4) then 2388 call hetransdipdens(2) !Finally call calcdelta_r 2389 else if (isel==5) then 2390 call exctransdip 2391 else if (isel==6) then 2392 call hetransdipdens(3) 2393 else if (isel==7) then 2394 write(*,"(a)") " Note: To calculate the ghost-hunter index proposed in JCC, 38, 2151 (2017), you should use option 1 of function 1 & 2395 to calculate hole-electron distribution, then the index will be automatically printed, followed by its two terms. See Section 3.21.7 of the manual for more details" 2396 write(*,*) "Press ENTER to contineu" 2397 read(*,*) 2398! write(*,"(a)") " PS: The index calculated in this way is somewhat different to the original paper, & 2399! in which the 1/D_CT term is calculated based on expensive relaxed density. If you really want to reproduce it, you can use function 3 & 2400! to calculate the 1/D_CT term corresponding to relaxed density, and then manually calculate ghost-hunter index" 2401 else if (isel==8) then 2402 call hetransdipdens(4) 2403 end if 2404 end do 2405 2406 2407!!!--------------------------------------- 2408!100!!------------------- Misc and some not important functions, Part 1 2409else if (infuncsel1==100) then 2410 do while(.true.) 2411 write(*,*) " ------------ Other functions (Part 1) ------------ " 2412 write(*,*) "0 Return" 2413 write(*,*) "1 Draw scatter graph between two functions and generate their cube files" 2414 write(*,*) "2 Export .pdb/.xyz/.wfn/.wfx/.molden/.fch/.47 files or input file of QC codes" 2415 write(*,*) "3 Calculate molecular van der Waals Volume" 2416 write(*,*) "4 Integrate a function in whole space" 2417 write(*,*) "5 Show overlap integral between alpha and beta orbitals" 2418 write(*,*) "6 Monitor SCF convergence process of Gaussian" 2419 write(*,*) "7 Generate Gaussian input file with initial guess from converged wavefunction" 2420 write(*,*) "8 Generate Gaussian input file with initial guess from fragment wavefunctions" 2421 write(*,*) "9 Evaluate coordination number for all atoms" 2422 write(*,*) "11 Calculate overlap and centroid distance between two orbitals" 2423 write(*,*) "13 Calculate HOMA and Bird aromaticity index" 2424 write(*,*) "14 Calculate LOLIPOP (LOL Integrated Pi Over Plane)" 2425 write(*,*) "15 Calculate intermolecular orbital overlap" 2426! write(*,*) "16 Calculate intermolecular tranfer integral by direct coupling method" 2427 write(*,*) "18 Yoshizawa's electron transport route analysis" 2428 write(*,*) "19 Generate promolecular .wfn file from fragment wavefunctions" 2429 write(*,*) "20 Calculate Hellmann-Feynman forces" 2430 write(*,*) "21 Calculate properties based on geometry information for specific atoms" 2431 write(*,*) "22 Detect pi orbitals and set occupation numbers" 2432 write(*,*) "23 Fit function distribution to atomic value" 2433 write(*,*) "24 Obtain NICS_ZZ for non-planar system" 2434 write(*,*) "25 Calculate area and perimeter for a ring" 2435 2436 read(*,*) infuncsel2 2437 if (infuncsel2==0) then 2438 goto 10 2439 else if (infuncsel2==1) then 2440 call funcvsfunc 2441 else if (infuncsel2==2) then 2442 write(*,*) "1 Output current structure to .pdb file" 2443 write(*,*) "2 Output current structure to .xyz file" 2444 write(*,*) "4 Output current wavefunction as .wfx file" 2445 write(*,*) "5 Output current wavefunction as .wfn file" 2446 write(*,*) "6 Output current wavefunction as Molden input file (.molden)" 2447 write(*,*) "7 Output current wavefunction as .fch file" 2448 write(*,*) "8 Output current wavefunction as .47 file" 2449 write(*,*) "10 Output current structure to Gaussian input file" 2450 write(*,*) "11 Output current structure to GAMESS-US input file" 2451 write(*,*) "12 Output current structure to ORCA input file" 2452 write(*,*) "13 Output current structure to NWChem input file" 2453 write(*,*) "14 Output current structure to MOPAC input file" 2454 write(*,*) "15 Output current structure to PSI input file" 2455 write(*,*) "16 Output current structure to MRCC input file" 2456 write(*,*) "17 Output current structure to CFOUR input file" 2457 write(*,*) "18 Output current structure to Molpro input file" 2458 write(*,*) "19 Output current structure to Dalton input file" 2459 write(*,*) "20 Output current structure to Molcas input file" 2460 read(*,*) itmp 2461 if (itmp==1) then 2462 write(*,*) "Input the path for pdb file, e.g. c:\ltwd.pdb" 2463 read(*,"(a)") c200tmp 2464 call outpdb(c200tmp,10) 2465 else if (itmp==2) then 2466 write(*,*) "Input the path for xyz file, e.g. c:\ltwd.xyz" 2467 read(*,"(a)") c200tmp 2468 call outxyz(c200tmp,10) 2469 else if (itmp==4) then 2470 if (.not.allocated(b)) then 2471 write(*,*) "Error: The input file you used does not contain GTF information!" 2472 else 2473 write(*,*) "Input the path, e.g. c:\ltwd.wfx" 2474 read(*,"(a)") c200tmp 2475 call outwfx(c200tmp,1,10) 2476 write(*,*) "Done!" 2477 end if 2478 else if (itmp==5) then 2479 if (.not.allocated(b)) then 2480 write(*,*) "Error: The input file you used does not contain GTF information!" 2481 else 2482 write(*,*) "Input the path, e.g. c:\ltwd.wfn" 2483 read(*,"(a)") c200tmp 2484 call outwfn(c200tmp,1,1,10) 2485 write(*,*) "Done!" 2486 end if 2487 else if (itmp==6) then 2488 if (.not.allocated(CObasa)) then 2489 write(*,*) "Error: This function works only when input file contains basis function information" 2490 else 2491 write(*,*) "Input the path, e.g. c:\ltwd.molden" 2492 read(*,"(a)") c200tmp 2493 write(*,*) "Exporting, please wait..." 2494 call outmolden(c200tmp,10) 2495 end if 2496 else if (itmp==7) then 2497 if (.not.allocated(CObasa)) then 2498 write(*,*) "Error: This function works only when input file contains basis function information" 2499 else 2500 write(*,*) "Input the path, e.g. c:\ltwd.fch" 2501 read(*,"(a)") c200tmp 2502 write(*,*) "Exporting, please wait..." 2503 call outfch(c200tmp,10) 2504 end if 2505 else if (itmp==8) then 2506 if (.not.allocated(CObasa)) then 2507 write(*,*) "Error: This function works only when input file contains basis function information" 2508 else 2509 write(*,*) "Input the path, e.g. c:\ltwd.47" 2510 read(*,"(a)") c200tmp 2511 write(*,*) "Exporting, please wait..." 2512 call out47(c200tmp,10) 2513 end if 2514 else if (itmp==10) then 2515 write(*,*) "Input the path, e.g. c:\ltwd.gjf" 2516 read(*,"(a)") c200tmp 2517 call outgjf(c200tmp,10) 2518 else if (itmp==11) then 2519 write(*,*) "Input the path, e.g. c:\ltwd.inp" 2520 read(*,"(a)") c200tmp 2521 call outGAMESSinp(c200tmp,10) 2522 else if (itmp==12) then 2523 write(*,*) "Input the path, e.g. c:\ltwd.inp" 2524 read(*,"(a)") c200tmp 2525 call outORCAinp(c200tmp,10) 2526 else if (itmp==13) then 2527 write(*,*) "Input the path, e.g. c:\ltwd.nw" 2528 read(*,"(a)") c200tmp 2529 call outNWCheminp(c200tmp,10) 2530 else if (itmp==14) then 2531 write(*,*) "Input the path, e.g. c:\ltwd.mop" 2532 read(*,"(a)") c200tmp 2533 call outMOPACinp(c200tmp,10) 2534 else if (itmp==15) then 2535 write(*,*) "Input the path, e.g. c:\ltwd.inp" 2536 read(*,"(a)") c200tmp 2537 call outPSIinp(c200tmp,10) 2538 else if (itmp==16) then 2539 c200tmp="MINP" 2540 call outMRCCinp(c200tmp,10) 2541 else if (itmp==17) then 2542 c200tmp="ZMAT" 2543 call outCFOURinp(c200tmp,10) 2544 else if (itmp==18) then 2545 write(*,*) "Input the path, e.g. c:\ltwd.inp" 2546 read(*,"(a)") c200tmp 2547 call outmolproinp(c200tmp,10) 2548 else if (itmp==19) then 2549 c200tmp=" " 2550 write(*,"(a)") " Input path of .dal file, e.g. c:\DFT.dal (directly press ENTER if you don't need it)" 2551 read(*,"(a)") c200tmp 2552 write(*,*) "Input path of .mol file, e.g. c:\ltwd.mol" 2553 read(*,"(a)") c200tmp2 2554 call outDaltoninp(c200tmp,c200tmp2,10) 2555 else if (itmp==20) then 2556 write(*,*) "Input the path, e.g. c:\ltwd.inp" 2557 read(*,"(a)") c200tmp 2558 call outmolcasinp(c200tmp,10) 2559 end if 2560 else if (infuncsel2==3) then 2561 if (MCvolmethod==1) then 2562 write(*,*) "100*2^i points will be used to evaluate the volume by Monte Carlo method" 2563 write(*,*) "The volume is defined as superposition of vdW sphere of atoms" 2564 write(*,*) 2565 do while(.true.) 2566 write(*,*) "Please input i, generally 10 is recommended, big system requires bigger i" 2567 write(*,*) "Input 0 can return" 2568 read(*,*) pointexp 2569 if (pointexp==0) exit 2570 call calcvolume(1,pointexp,0D0,1D0) 2571 end do 2572 else if (MCvolmethod==2) then 2573 write(*,*) "100*2^i points will be used to evaluate the volume by Monte Carlo method." 2574 write(*,"(a)") " The volume is defined as the region encompassed by the isosurface of density equals to x, & 2575 The box used in Monte Carlo procedure will be enlarged by k multiples vdW radius in each side." 2576 write(*,"(a)") " Hint: For evaluating the volume encompassed by 0.001 isosurface of density for small molecule, & 2577 we suggest you input 9,0.001,1.7" 2578 write(*,*) 2579 do while(.true.) 2580 write(*,*) "Please input i,x,k (Input 0,0,0 can return)" 2581 read(*,*) pointexp,tmpisoval,enlarbox 2582 if (pointexp==0.and.tmpisoval==0.and.enlarbox==0) exit 2583 call calcvolume(2,pointexp,tmpisoval,enlarbox) 2584 end do 2585 end if 2586 else if (infuncsel2==4) then 2587 if (ispecial/=1) then 2588 call selfunc_interface(ifunc) 2589 call intfunc(ifunc) 2590 else if (ispecial==1) then 2591 call intfunc(1) 2592 end if 2593 else if (infuncsel2==-4) then 2594 call intdiff(1) 2595 else if (infuncsel2==-5) then 2596 call intdiff(2) 2597 else if (infuncsel2==5) then 2598 call aboverlap 2599 else if (infuncsel2==6) then 2600 call monitorscf 2601 else if (infuncsel2==7) then 2602 call gengauguess 2603 else if (infuncsel2==8) then 2604 call fragguess 2605 else if (infuncsel2==9) then 2606 call coordnum 2607 else if (infuncsel2==11) then 2608 call ovlpdistorb 2609 else if (infuncsel2==13) then 2610 call HOMA_Bird 2611 else if (infuncsel2==14) then 2612 call LOLIPOP 2613 else if (infuncsel2==15) then 2614 call intmolovlp 2615 else if (infuncsel2==16) then 2616 call intmoltransint 2617 else if (infuncsel2==18) then 2618 call Yoshieletrans 2619 else if (infuncsel2==19) then 2620 call genpromolwfn 2621 else if (infuncsel2==20) then 2622 call hellmann_feynman 2623 else if (infuncsel2==21) then 2624 call calcgeomprop 2625 else if (infuncsel2==22) then 2626 call procpiorb 2627 else if (infuncsel2==23) then 2628 call fitfunc 2629 else if (infuncsel2==24) then 2630 call utilNICS_ZZ 2631 else if (infuncsel2==25) then 2632 call calcringsize 2633 end if 2634 write(*,*) 2635 end do 2636 2637 2638!!!--------------------------------------- 2639!200!!------------------- Misc and not important functions, Part 2 2640else if (infuncsel1==200) then 2641 do while(.true.) 2642 write(*,*) " ------------ Other functions (Part 2) ------------ " 2643 write(*,*) "0 Return" 2644 write(*,*) "1 Weak interaction analysis for fluctuation environment by RDG method" 2645 write(*,*) "2 Calculate atomic and bond dipole moments in Hilbert space" 2646 write(*,*) "3 Generate cube file for multiple orbital wavefunctions" 2647 write(*,*) "4 Generate iso-chemical shielding surfaces (ICSS) and related quantities" 2648 write(*,*) "5 Plot radial distribution function for a real space function" 2649 write(*,*) "6 Analyze correspondence between orbitals in two wavefunctions" 2650 write(*,*) "7 Parse output of (hyper)polarizability task of Gaussian" 2651 write(*,*) "8 Calculate (hyper)polarizability by sum-over-states (SOS) method" 2652 write(*,*) "9 Calculate average bond length and average coordinate number" 2653 write(*,*) "10 Output various kinds of integral between orbitals" 2654 write(*,*) "11 Calculate center, the first and second moments of a real space function" 2655 write(*,*) "12 Calculate energy index (EI) or bond polarity index (BPI)" 2656 write(*,*) "13 Pipek-Mezey orbital localization" 2657 write(*,*) "14 Perform integration within isosurfaces of a real space function" 2658 write(*,*) "15 Calculate electron correlation index (PCCP, 18, 24015)" 2659 write(*,*) "16 Generate natural orbitals based on the density matrix in .fch/.fchk file" 2660! write(*,*) "20 Calculate electronic coupling matrix element by FCD or GMH method" 2661 read(*,*) infuncsel2 2662 if (infuncsel2==0) then 2663 goto 10 2664 else if (infuncsel2==1) then 2665 call RDG_MD 2666 else if (infuncsel2==2) then 2667 call atmbonddip 2668 else if (infuncsel2==3) then 2669 call genmultiorbcube 2670 else if (infuncsel2==4) then 2671 call ICSS 2672 else if (infuncsel2==5) then 2673 call plotraddis 2674 else if (infuncsel2==6) then 2675 call orbcorres 2676 else if (infuncsel2==7) then 2677 call parseGauPolar 2678 else if (infuncsel2==8) then 2679 call SOS 2680 else if (infuncsel2==9) then 2681 call atmavgdist 2682 else if (infuncsel2==10) then 2683 call outorbint 2684 else if (infuncsel2==11) then 2685 call funcmoment 2686 else if (infuncsel2==12) then 2687 call calcEIBPI 2688 else if (infuncsel2==13) then 2689 call pipek_mezey 2690 else if (infuncsel2==14) then 2691 call intisosurface 2692 else if (infuncsel2==15) then 2693 call elecorridx 2694 else if (infuncsel2==16) then 2695 call fch_gennatorb 2696! else if (infuncsel2==20) then 2697! call FCD 2698 end if 2699 write(*,*) 2700 end do 2701 2702 2703else if (infuncsel1==98) then 2704 call sphatmraddens 2705 2706 2707!!!--------------------------------------- 2708!1000!!------------------- Special functions 2709else if (infuncsel1==1000) then 27101000 write(*,*) 2711 write(*,*) "0 Return to main menu" 2712 write(*,"(a,3f12.6)") " 1 Set reference point, current(Bohr):",refx,refy,refz 2713 write(*,"(a,i6)") " 2 Set iuserfunc, current:",iuserfunc 2714 write(*,"(a,i6)") " 3 Set iskipnuc, current:",iskipnuc 2715 if (pleA==0D0.and.pleB==0D0.and.pleC==0D0.and.pleD==0D0) then 2716 write(*,"(a)") " 4 Set the plane for user-defined function 38 (Not defined)" 2717 else 2718 write(*,"(a)") " 4 Set the plane for user-defined function 38 (Defined)" 2719 end if 2720 write(*,"(a,1PD18.8)") " 5 Set global temporary variable, current:",globaltmp 2721 write(*,"(a,i3)") " 10 Set the number of threads, current:", getNThreads() 2722 write(*,*) "90 Calculate nuclear attractive energy between a fragment and an orbital" 2723 write(*,*) "97 Generate natural orbitals based on density matrix outputted by MRCC program" 2724 write(*,*) "98 Generate natural orbitals based on density matrix in .fch/.fchk" 2725 write(*,*) "99 Show EDF information (if any)" 2726 write(*,*) "100 Check the sanity of present wavefunction" 2727 read(*,*) i 2728 if (i==1) then 2729 write(*,*) "Input x,y,z in Bohr, e.g. 3.0,0.0,1.3" 2730 read(*,*) refx,refy,refz 2731 write(*,*) "Done!" 2732 else if (i==2) then 2733 write(*,*) "Input an integer, e.g. 24" 2734 read(*,*) iuserfunc 2735 write(*,*) "Done!" 2736 else if (i==3) then 2737 write(*,*) "Input the index of the nucleus, e.g. 24" 2738 read(*,*) iskipnuc 2739 write(*,*) "Done!" 2740 else if (i==4) then 2741 write(*,*) "1 Input index of three atoms to define the plane" 2742 write(*,*) "2 Input XYZ coordinate of three points to define the plane" 2743 read(*,*) iseldef 2744 if (iseldef==1) then 2745 write(*,*) "Input three indices, e.g. 2,4,5" 2746 read(*,*) i1,i2,i3 2747 call pointABCD(a(i1)%x,a(i1)%y,a(i1)%z,a(i2)%x,a(i2)%y,a(i2)%z,a(i3)%x,a(i3)%y,a(i3)%z,pleA,pleB,pleC,pleD) 2748 else if (iseldef==2) then 2749 write(*,*) "Input coordinate for point 1 (in Bohr), e.g. 1.0,-0.2,0.3" 2750 read(*,*) xtmp1,ytmp1,ztmp1 2751 write(*,*) "Input coordinate for point 2 (in Bohr), e.g. 2.0,-0.3,0.1" 2752 read(*,*) xtmp2,ytmp2,ztmp2 2753 write(*,*) "Input coordinate for point 3 (in Bohr), e.g. 1.3,-1.2,0.33" 2754 read(*,*) xtmp3,ytmp3,ztmp3 2755 call pointABCD(xtmp1,ytmp1,ztmp1,xtmp2,ytmp2,ztmp2,xtmp3,ytmp3,ztmp3,pleA,pleB,pleC,pleD) 2756 end if 2757 tmpval=dsqrt(pleA**2+pleB**2+pleC**2) 2758 write(*,"(' The unit vector normal to the plane is:',3f10.5)") pleA/tmpval,pleB/tmpval,pleC/tmpval 2759 goto 1000 2760 else if (i==5) then 2761 write(*,*) "Input the value" 2762 read(*,*) globaltmp 2763 else if (i==10) then 2764 write(*,*) "Input an integer, e.g. 8" 2765 read(*,*) iniNThreads 2766 write(*,*) "Done!" 2767 else if (i==90) then 2768 call attene_orb_fragnuc 2769 else if (i==97) then 2770 call MRCC_gennatorb 2771 else if (i==98) then 2772 call fch_gennatorb 2773 else if (i==99) then 2774 if (.not.allocated(b_EDF)) then 2775 write(*,*) "EDF field was not loaded" 2776 else 2777 write(*,"( ' The number of inner-core electrons represented by EDF:',i6)") nEDFelec 2778 write(*,*) "Information of EDF primitives:" 2779 write(*,*) "Column 1: Index" 2780 write(*,*) "Column 2: Atom" 2781 write(*,*) "Column 3: Function type" 2782 write(*,*) "Column 4: Exponent" 2783 write(*,*) "Column 5: Coefficient" 2784 do iEDFprim=1,nEDFprims 2785 write(*,"(3i6,2f20.8)") iEDFprim,b_EDF(iEDFprim)%center,b_EDF(iEDFprim)%functype,b_EDF(iEDFprim)%exp,CO_EDF(iEDFprim) 2786 end do 2787 end if 2788 else if (i==100) then 2789 call wfnsanity 2790 end if 2791 2792end if 2793 2794end do !End main cycle 2795 2796 2797 2798 2799contains 2800!=================================================================================================== 2801!=================================================================================================== 2802!!!!!!!!!!!!!!!! Below routines are contained in and closely related to main program !!!!!!!!!!!!!!! 2803!=================================================================================================== 2804!=================================================================================================== 2805 2806 2807 2808 2809!!!-------------------------- Set content of custom plot 2810subroutine customplotsetup 2811integer i 2812write(*,*) "How many files to deal with? (Excluding the file that has been loaded)" 2813read(*,*) ncustommap 2814if (allocated(custommapname)) deallocate(custommapname) 2815if (allocated(customop)) deallocate(customop) 2816allocate(custommapname(ncustommap)) 2817allocate(customop(ncustommap)) 2818write(*,*) "The avaliable operator: +,-,*,/" !* can also be intputted as x 2819write(*,*) "e.g. -,sob.wfn means minus the property of sob.wfn from the first file" 2820do i=1,ncustommap 2821 write(*,"('Input the operator and filename of',i5)") i 2822 do while(.true.) 2823 read(*,"(a1,1x,a)") customop(i),custommapname(i) 2824 inquire(file=custommapname(i),exist=alive) 2825 if (alive.eqv..true.) then 2826 exit 2827 else 2828 write(*,*) "File not found, input again" 2829 end if 2830 end do 2831end do 2832end subroutine 2833 2834 2835!!!----------------------Set contour line 2836subroutine setctr 2837implicit real*8 (a-h,o-z) 2838character outfilename*80,selectyn*1 2839do while(.true.) 2840 if (j/=6) then !If last operation is not saving file 2841 write(*,*) "Current isovalue line:" 2842 do j=1,lastctrval !The lastest contour line serial 2843 write(*,"(i3,':',f19.8,' ')",advance='no') j,ctrval(j) 2844 if (mod(j,3)==0) write(*,*) 2845 end do 2846 if (mod(lastctrval,3)/=0) write(*,*) 2847 end if 2848 if (lastctrval==0) write(*,*) "None" 2849 if (allocated(boldlinelist)) write(*,*) "Bolded line index:",boldlinelist 2850 write(*,*) "1 Save setting and return" 2851 write(*,*) "2 Replace a value of contour line by inputting" 2852 write(*,*) "3 Add a new contour line" 2853 write(*,*) "4 Delete a contour line" 2854 write(*,*) "5 Delete a range of contour lines" 2855 write(*,*) "6 Save contour setting to external file" 2856 write(*,*) "7 Load contour setting from external file" 2857 write(*,*) "8 Generate contour value by arithmetic progression" 2858 write(*,*) "9 Generate contour value by geometric series" 2859 if (.not.allocated(boldlinelist)) write(*,*) "10 Enable some contour lines bold style" 2860 if (allocated(boldlinelist)) write(*,*) "10 Disable some contour lines bold style" 2861 write(*,*) "11 Set color for positive contour lines" 2862 write(*,*) "12 Set line style and width for positive contour lines" 2863 write(*,*) "13 Set color for negative contour lines" 2864 write(*,*) "14 Set line style and width for negative contour lines" 2865 write(*,*) "15 Set line style and width suitable for publication" 2866 read(*,*) j 2867 if (j==1) then 2868 exit 2869 else if (j==2) then 2870 write(*,*) "Replace which value of contour line by what value?" 2871 write(*,*) "(e.g. Input 4,0.015 to replace the fourth value of contour line by 0.015)" 2872 read(*,*) j,selfctrval 2873 if (j>=1.and.j<=lastctrval) then 2874 ctrval(j)=selfctrval 2875 else 2876 write(*,*) "The number exceed valid range" 2877 end if 2878 else if (j==3) then 2879 if (lastctrval<size(ctrval)) then 2880 write(*,*) "Input the value of contour line you want to add" 2881 read(*,*) selfctrval 2882 lastctrval=lastctrval+1 2883 ctrval(lastctrval)=selfctrval 2884 else 2885 write(*,"(a,i8)") " Error: The number of contour line couldn't exceed",size(ctrval) 2886 end if 2887 else if (j==4) then 2888 write(*,*) "Delete which contour line? Input a number" 2889 read(*,*) j 2890 if (j>=1.and.j<=lastctrval) then 2891 lastctrval=lastctrval-1 2892 ctrval(j:lastctrval)=ctrval(j+1:lastctrval+1) 2893 else 2894 write(*,*) "The number exceed valid range" 2895 end if 2896 else if (j==5) then 2897 write(*,*) "Delete the range of contour line number (e.g. 3,10)" 2898 write(*,*) "Note: You can also input 0,0 to delete all" 2899 read(*,*) j1,j2 2900 if (j1==0.and.j2==0) then 2901 lastctrval=0 2902 else if (j2>j1.and.j1>=1.and.j2<=lastctrval) then 2903 lastctrval=lastctrval-(j2-j1+1) 2904 ctrval(j1:lastctrval)=ctrval(j2+1:lastctrval+(j2-j1+1)) 2905 else 2906 write(*,*) "Invalid input, input excced current range" 2907 end if 2908 else if (j==6) then 2909 write(*,*) "Input the filename for outputting current setting e.g. c:\ltwd.txt" 2910 read(*,*) outfilename 2911 open(10,file=outfilename) 2912 write(10,"(f19.8)") (ctrval(i),i=1,lastctrval) 2913 close(10) 2914 write(*,"(a,a)") " Contour setting has been saved to ",trim(outfilename) 2915 write(*,*) 2916 else if (j==7) then 2917 write(*,*) "Input filename, e.g. C:\ctr.txt" 2918 do while(.true.) 2919 read(*,"(a)") extctrsetting 2920 inquire(file=extctrsetting,exist=alive) 2921 if (alive .eqv. .true.) exit 2922 write(*,*) "File not found, input again" 2923 end do 2924 open(10,file=extctrsetting,access="sequential",status="old") 2925 ierror=0 2926 lastctrval=0 2927 do i=1,size(ctrval) 2928 read(10,*,iostat=ierror) temp 2929 if (ierror/=0) exit 2930 ctrval(i)=temp 2931 lastctrval=i 2932 end do 2933 close(10) 2934 else if (j==8.or.j==9) then 2935 write(*,*) "Input start value, step and total number" 2936 if (j==8) write(*,*) "e.g. 0.4,0.1,10 generate 0.4,0.5,0.6,0.7 ... 1.3" 2937 if (j==9) write(*,*) "e.g. 2,3,10 generate 2,6,18,54 ... 39366" 2938 read(*,*) ctrgenstrval,ctrgenstep,igennum 2939 write(*,*) "If remove existed contour lines? (y/n)" 2940 read(*,*) selectyn 2941 if (selectyn=='y'.or.selectyn=='Y') then 2942 lastctrval=0 2943 else if (selectyn=='n'.or.selectyn=='N') then !Append to existed contour lines 2944 if (lastctrval+igennum>size(ctrval)) then 2945 igennum=size(ctrval)-lastctrval 2946 write(*,*) "Warning: The total number of coutour lines exceeded the upper limit!" 2947 end if 2948 end if 2949 do jj=1,igennum 2950 if (j==8) ctrval(lastctrval+jj)=ctrgenstrval+ctrgenstep*(jj-1) 2951 if (j==9) ctrval(lastctrval+jj)=ctrgenstrval*ctrgenstep**(jj-1) 2952 end do 2953 lastctrval=lastctrval+igennum 2954 else if (j==10) then 2955 if (allocated(boldlinelist)) then !Disable bold line 2956 deallocate(boldlinelist) 2957 write(*,*) "No line is bolded now, you can select this function again to set bolded lines" 2958 else 2959 write(*,*) "Set how many lines bold?" 2960 read(*,*) numboldline 2961 allocate(boldlinelist(numboldline)) 2962 do itmp=1,numboldline 2963 write(*,"(' Input the index of bolded line',i4)") itmp 2964 read(*,*) boldlinelist(itmp) 2965 end do 2966 end if 2967 else if (j==11) then 2968 write(*,*) "Use which color?" 2969 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 2970 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 2971 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 2972 read(*,*) iclrindctrpos 2973 else if (j==12) then 2974 write(*,*) "Input length of line segment and interstice" 2975 write(*,*) "e.g. 1,0 means solid line; 1,10 means DOT; 10,10 means DASH" 2976 write(*,*) " 10,15 means DASH with larger interstice" 2977 write(*,*) "Note: 1,0 and 10,15 are default for positive and negative lines, respectively" 2978 read(*,*) ctrposstyle(1),ctrposstyle(2) 2979 write(*,*) "Input line width, e.g. 2" 2980 read(*,*) iwidthposctr 2981 else if (j==13) then 2982 write(*,*) "Use which color?" 2983 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 2984 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 2985 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 2986 read(*,*) iclrindctrneg 2987 else if (j==14) then 2988 write(*,*) "Input length of line segment and interstice" 2989 write(*,*) "e.g. 1,0 means solid line; 1,10 means DOT; 10,10 means DASH" 2990 write(*,*) " 10,15 means DASH with larger interstice" 2991 write(*,*) "Note: 1,0 and 10,15 are default for positive and negative lines, respectively" 2992 read(*,*) ctrnegstyle(1),ctrnegstyle(2) 2993 write(*,*) "Input line width, e.g. 2" 2994 read(*,*) iwidthnegctr 2995 else if (j==15) then 2996 iclrindctrpos=11 2997 iclrindctrneg=3 2998 ctrposstyle(1)=1 2999 ctrposstyle(2)=0 3000 iwidthposctr=6 3001 ctrnegstyle(1)=10 3002 ctrnegstyle(2)=15 3003 iwidthnegctr=6 3004 write(*,"(a)") " Done. The saved picture with current line setting is suitable for publication purpose" 3005 end if 3006end do 3007end subroutine 3008 3009 3010!!------------------ Set marker of CPs and paths on contour/gradient map 3011subroutine settopomark 3012implicit real*8 (a-h,o-z) 3013do while(.true.) 3014 write(*,*) "0 Return" 3015 if (imark3n3==0) write(*,*) "1 Show (3,-3) CPs" 3016 if (imark3n3==1) write(*,*) "1 Don't show (3,-3) CPs" 3017 if (imark3n1==0) write(*,*) "2 Show (3,-1) CPs" 3018 if (imark3n1==1) write(*,*) "2 Don't show (3,-1) CPs" 3019 if (imark3p1==0) write(*,*) "3 Show (3,+1) CPs" 3020 if (imark3p1==1) write(*,*) "3 Don't show (3,+1) CPs" 3021 if (imark3p3==0) write(*,*) "4 Show (3,+3) CPs" 3022 if (imark3p3==1) write(*,*) "4 Don't show (3,+3) CPs" 3023 if (imarkpath==0) write(*,*) "5 Show paths" 3024 if (imarkpath==1) write(*,*) "5 Don't show paths" 3025 write(*,*) "10 Set size of markers of CPs" 3026 write(*,*) "11 Set thickness of paths" 3027 write(*,*) "12 Set color of paths" 3028 write(*,*) "13 Set thickness of the interbasin paths derived from (3,-1)" 3029 write(*,*) "14 Set color of the interbasin paths derived from (3,-1)" 3030 read(*,*) isel 3031 3032 if (isel==0) then 3033 exit 3034 else if (isel==1) then 3035 if (imark3n3==1) then 3036 imark3n3=0 3037 else 3038 imark3n3=1 3039 end if 3040 else if (isel==2) then 3041 if (imark3n1==1) then 3042 imark3n1=0 3043 else 3044 imark3n1=1 3045 end if 3046 else if (isel==3) then 3047 if (imark3p1==1) then 3048 imark3p1=0 3049 else 3050 imark3p1=1 3051 end if 3052 else if (isel==4) then 3053 if (imark3p3==1) then 3054 imark3p3=0 3055 else 3056 imark3p3=1 3057 end if 3058 else if (isel==5) then 3059 if (imarkpath==1) then 3060 imarkpath=0 3061 else 3062 imarkpath=1 3063 end if 3064 else if (isel==10) then 3065 write(*,*) "Input a value, e.g. 30" 3066 read(*,*) sizemarkcp 3067 else if (isel==11) then 3068 write(*,*) "Input a value, e.g. 5" 3069 read(*,*) sizemarkpath 3070 else if (isel==12) then 3071 write(*,*) "Input R,G,B value, between 0.0 to 1.0. e.g. 0.5,0.7.1.0" 3072 write(*,*) "Note: 0,0,0 corresponds to black, 1,1,1 corresponds to white" 3073 read(*,*) clrRpath,clrGpath,clrBpath 3074 else if (isel==13) then 3075 write(*,*) "Input a value, e.g. 5" 3076 read(*,*) sizemark3n1path 3077 else if (isel==14) then 3078 write(*,*) "Input R,G,B value, between 0.0 to 1.0. e.g. 0.5,0.7.1.0" 3079 write(*,*) "Note: 0,0,0 corresponds to black, 1,1,1 corresponds to white" 3080 read(*,*) clrR3n1path,clrG3n1path,clrB3n1path 3081 end if 3082end do 3083end subroutine 3084 3085 3086 3087 3088end program 3089