1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! Copyright 2010. Los Alamos National Security, LLC. This material was ! 3! produced under U.S. Government contract DE-AC52-06NA25396 for Los Alamos ! 4! National Laboratory (LANL), which is operated by Los Alamos National ! 5! Security, LLC for the U.S. Department of Energy. The U.S. Government has ! 6! rights to use, reproduce, and distribute this software. NEITHER THE ! 7! GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, ! 8! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS ! 9! SOFTWARE. If software is modified to produce derivative works, such ! 10! modified software should be clearly marked, so as not to confuse it ! 11! with the version available from LANL. ! 12! ! 13! Additionally, this program is free software; you can redistribute it ! 14! and/or modify it under the terms of the GNU General Public License as ! 15! published by the Free Software Foundation; version 2.0 of the License. ! 16! Accordingly, this program is distributed in the hope that it will be ! 17! useful, but WITHOUT ANY WARRANTY; without even the implied warranty of ! 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General ! 19! Public License for more details. ! 20!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 21 22MODULE LATTE_LIB 23 24 USE CONSTANTS_MOD 25 USE TIMER_MOD 26 USE SETUPARRAY 27 USE PPOTARRAY 28 USE PUREARRAY 29 USE COULOMBARRAY 30 USE SPINARRAY 31 USE SPARSEARRAY 32 USE MDARRAY 33 USE MYPRECISION 34 USE VIRIALARRAY 35 USE DIAGARRAY 36 USE KSPACEARRAY 37 USE LATTEPARSER_LATTE_MOD 38 USE NEBLISTARRAY 39 USE NONOARRAY 40 USE CONSTRAINTS_MOD 41 42#ifdef PROGRESSON 43 USE PRG_SYSTEM_MOD ! FROM PROGRESS 44 USE PRG_PULAYMIXER_MOD 45 USE PRG_EXTRAS_MOD 46 USE MIXER_MOD 47 USE BML 48#endif 49 50#ifdef MPI_ON 51 USE MPI 52#endif 53#ifdef DBCSR_ON 54 USE DBCSR_VAR_MOD 55#endif 56 57 IMPLICIT NONE 58 59 PRIVATE 60 61 ! Defines the version of the binary interface. 62 ! Adjust if non-backward compatible changes are made to the interface. 63 ! Use in codes using the library interface to make certain a compatible 64 ! version of the LATTE library is used. 65 INTEGER, PARAMETER :: LATTE_ABIVERSION = 20180622 66 67 PUBLIC :: LATTE, LATTE_ABIVERSION 68 69CONTAINS 70 71 !! Latte subroutine 72 !! \param FLAGS Different control flags that can be passed to LATTE (not in use yet) 73 !! \param NATS Number of atoms 74 !! \param COORDS Coordinates. Example: y-coordinate of atom 1 = COORDS(2,1) 75 !! \param TYPES An index for all the different atoms in the system. 76 !! \param NTYPES Number of different elements in the system 77 !! \param MASSES Element masses for every different element of the system. 78 !! \param XLO Lowest dimensions of the box 79 !! \param XHI Highest dimensions of the box 80 !! \param XY Tilt factor. 81 !! \param XZ Tilt factor. 82 !! \param YZ Tilt factor. The lattice vectors are constructed as: 83 !! a = (xhi-xlo,0,0); b = (xy,yhi-ylo,0); c = (xz,yz,zhi-zlo). 84 !! \param FORCES Forces for every atom as output. 85 !! \param MAXITER Latte MAXITER keyword. If MAXITER = -1, only the Forces are computed. 86 !! If MAXITER = 0, MAXITER is read from latte.in file. 87 !! IF MAXITER > 0, MAXITER is passed trough the library call. 88 !! \param VENERG This is the potential Energy that is given back from latte to the hosting code. 89 !! \param VEL Velocities passed to latte. 90 !! \param DT integration step passed to latte. 91 !! \param VIRIAL_INOUT Components of the second virial coefficient. 92 !! \param NEWSYSTEM Tells LATTE if a new system is passed. 93 !! \param EXISTERROR_INOUT Returns an error flag (.true.) to the hosting code. 94 !! 95 !! \brief This routine will be used load call latte_lib from a C/C++ program: 96 !! 97 !! \brief Note: To get the mass of atom 3 we do: 98 !! \verbatim 99 !! MASS(TYPES(3)) 100 !! \endverbatim 101 !! 102 !! \brief Note: To get the lattice vectors as formated in LATTE we do: 103 !! \verbatim 104 !! BOX(1,1) = XHI(1) - XLO(1); ... 105 !! \endverbatim 106 !! 107 !! \brief Note: All units are LATTE units by default. See https://github.com/losalamos/LATTE/blob/master/Manual/LATTE_manual.pdf 108 !! 109 SUBROUTINE LATTE(NTYPES, TYPES, CR_IN, MASSES_IN, XLO, XHI, XY, XZ, YZ, FTOT_OUT, & 110 MAXITER_IN, VENERG, VEL_IN, DT_IN, VIRIAL_INOUT, NEWSYSTEM, EXISTERROR_INOUT) 111 112 USE CONSTANTS_MOD, ONLY: EXISTERROR 113 114 IMPLICIT NONE 115 116 INTEGER :: I, J 117 INTEGER :: START_CLOCK, STOP_CLOCK, CLOCK_RATE, CLOCK_MAX 118 INTEGER :: MYID = 0 119 REAL :: TARRAY(2), RESULT, SYSTDIAG, SYSTPURE 120 REAL(LATTEPREC) :: DBOX 121 CHARACTER(LEN=50) :: FLNM 122 123 REAL(LATTEPREC), INTENT(IN) :: CR_IN(:,:),VEL_IN(:,:), MASSES_IN(:),XLO(3),XHI(3) 124 REAL(LATTEPREC), INTENT(IN) :: DT_IN, XY, XZ, YZ 125 REAL(LATTEPREC), INTENT(OUT) :: FTOT_OUT(:,:), VENERG 126 REAL(LATTEPREC), INTENT(OUT) :: VIRIAL_INOUT(6) 127 INTEGER, INTENT(IN) :: NTYPES, TYPES(:), MAXITER_IN 128 LOGICAL(1), INTENT(INOUT) :: EXISTERROR_INOUT 129 REAL(LATTEPREC) :: MLSI, LUMO, HOMO 130 INTEGER, INTENT(INOUT) :: NEWSYSTEM 131 132#ifdef PROGRESSON 133 TYPE(SYSTEM_TYPE) :: SY 134#endif 135 136#ifdef MPI_ON 137 INTEGER :: IERR, STATUS(MPI_STATUS_SIZE), NUMPROCS 138 CALL MPI_INIT( IERR ) 139 CALL MPI_COMM_RANK( MPI_COMM_WORLD, MYID, IERR ) 140 CALL MPI_COMM_SIZE( MPI_COMM_WORLD, NUMPROCS, IERR ) 141#endif 142 143 EXISTERROR = .FALSE. !We assume we start the lib call without errors 144 145 IF(.NOT. LIBINIT .OR. NEWSYSTEM == 1)THEN 146 147 CALL DEALLOCATEALL() 148 149 LIBRUN = .TRUE. 150 151 LIBCALLS = 0 ; MAXITER = -10 152 153 ! Only LATTE main code will create the animate folder 154 ! INQUIRE( FILE="animate/.", exist=LATTEINEXISTS) 155 ! IF (.NOT. LATTEINEXISTS) CALL SYSTEM("mkdir animate") 156 157 NUMSCF = 0 158 CHEMPOT = ZERO 159 160 ! Start timers 161 TX = INIT_TIMER() 162 TX = START_TIMER(LATTE_TIMER) 163 164 INQUIRE( FILE=LATTEINNAME, exist=LATTEINEXISTS ) 165 166 IF (LATTEINEXISTS) THEN 167 CALL PARSE_CONTROL(LATTEINNAME) 168 169#ifdef PROGRESSON 170 CALL PRG_PARSE_MIXER(MX,LATTEINNAME) 171#endif 172 173 ELSE 174 CALL READCONTROLS 175 ENDIF 176 177 IF(VERBOSE <= 0)THEN 178 OPEN(UNIT=6, FILE="/dev/null", FORM="formatted") 179 ELSE 180 OPEN(UNIT=6, FILE="log.latte", FORM="formatted") 181 ENDIF 182 183 IF(VERBOSE >= 1)THEN 184 WRITE(*,*)"# The log file for latte_lib" 185 WRITE(*,*)"" 186 CALL TIMEDATE_TAG("LATTE started at : ") 187 ENDIF 188 189 CALL READTB 190 191 IF (RESTART .EQ. 0) THEN 192 193 BOX = 0.0d0 194 BOX(1,1) = xhi(1) - xlo(1) 195 BOX(2,1) = XY 196 BOX(2,2) = xhi(2) - xlo(2) 197 BOX(3,1) = XZ 198 BOX(3,2) = YZ 199 BOX(3,3) = xhi(3) - xlo(3) 200 201 IF(VERBOSE >= 1)THEN 202 WRITE(*,*)"Lattice vectors:" 203 WRITE(*,*)"a=",BOX(1,1),BOX(1,2),BOX(1,3) 204 WRITE(*,*)"b=",BOX(2,1),BOX(2,2),BOX(2,3) 205 WRITE(*,*)"c=",BOX(3,1),BOX(3,2),BOX(3,3) 206 WRITE(*,*)"" 207 ENDIF 208 209 BOX_OLD = BOX 210 211 NATS = SIZE(CR_IN,DIM=2) 212 213 IF (.NOT.ALLOCATED(CR)) ALLOCATE(CR(3,NATS)) 214 CR = CR_IN 215 216 IF(VERBOSE >= 1)WRITE(*,*)"Converting masses to symbols ..." 217 IF(.NOT. ALLOCATED(ATELE)) ALLOCATE(ATELE(NATS)) 218 CALL MASSES2SYMBOLS(TYPES,NTYPES,MASSES_IN,NATS,ATELE) 219 220 !Forces, charges and element pointers are allocated in readcr 221 CALL READCR 222 223 FLUSH(6) 224 225 ELSE 226 227 IF(VERBOSE >= 1)WRITE(*,*)"Restarting calculation from file ..." 228 CALL READRESTART 229 230 ENDIF 231 232 IF (VERBOSE >= 1) WRITE(*,*)"Reading ppots from file (if PPOTON >= 1) ..." 233 IF (PPOTON .EQ. 1) CALL READPPOT 234 IF (PPOTON .EQ. 2) CALL READPPOTTAB 235 IF (PPOTON .EQ. 3) CALL READPPOTSPLINE 236 237 IF (DEBUGON .EQ. 1) THEN 238 CALL PLOTUNIV 239 IF (PPOTON .EQ. 1) CALL PLOTPPOT 240 ENDIF 241 242 CALL GETHDIM 243 244 CALL GETMATINDLIST 245 246 IF (VERBOSE >= 1) WRITE(*,*)"Getting rho0 ..." 247 CALL RHOZERO 248 249 CALL GETBNDFIL() 250 251 FLUSH(6) 252 253#ifdef GPUON 254 255 CALL INITIALIZE( NGPU ) 256 257#endif 258 259#ifdef DBCSR_ON 260 261 IF (CONTROL .EQ. 2 .AND. SPARSEON .EQ. 1) CALL INIT_DBCSR 262 263#endif 264 265 IF (KBT .LT. 0.0000001 .OR. CONTROL .EQ. 2) ENTE = ZERO 266 267 IF (.NOT. ALLOCATED(V)) ALLOCATE(V(3,NATS)) 268 IF(VERBOSE >= 1)WRITE(*,*)"End of INITIALIZATION" 269 270 ELSE 271 272 BOX = 0.0d0 273 BOX(1,1) = xhi(1) - xlo(1) 274 BOX(2,1) = XY 275 BOX(2,2) = xhi(2) - xlo(2) 276 BOX(3,1) = XZ 277 BOX(3,2) = YZ 278 BOX(3,3) = xhi(3) - xlo(3) 279 280 IF(VERBOSE >= 1)THEN 281 WRITE(*,*)"Lattice vectors:" 282 WRITE(*,*)"a=",BOX(1,1),BOX(1,2),BOX(1,3) 283 WRITE(*,*)"b=",BOX(2,1),BOX(2,2),BOX(2,3) 284 WRITE(*,*)"c=",BOX(3,1),BOX(3,2),BOX(3,3) 285 WRITE(*,*)"" 286 ENDIF 287 288 LIBCALLS = LIBCALLS + 1 289 290 NATS = SIZE(CR_IN,DIM=2) 291 292 IF(.NOT.ALLOCATED(CR)) ALLOCATE(CR(3,NATS)) 293 CR = CR_IN 294 295 FLUSH(6) 296 297 ENDIF 298 !End of initialization 299 300 301 IF (MDON .EQ. 0 .AND. RELAXME .EQ. 0 .AND. DOSFITON .EQ. 0 & 302 .AND. PPFITON .EQ. 0 .AND. ALLFITON .EQ. 0) THEN 303 304 305 ! IF (.NOT. LIBINIT) THEN 306 307 ! 308 ! Start the timers 309 ! 310 311 CALL SYSTEM_CLOCK(START_CLOCK, CLOCK_RATE, CLOCK_MAX) 312 313#ifndef FCIDxlf 314 CALL DTIME(TARRAY, RESULT) 315#endif 316 317 318 ! Set up neighbor lists for building the H and pair potentials 319 320 CALL ALLOCATENEBARRAYS 321 322 323 IF (ELECTRO .EQ. 1) THEN 324 325 CALL ALLOCATECOULOMB 326 327 CALL INITCOULOMB 328 329 ENDIF 330 331 332 IF (BASISTYPE .EQ. "NONORTHO") CALL ALLOCATENONO 333 334 335 336 CALL NEBLISTS(0) 337 338 339 ! Build the charge independent H matrix 340 341 IF (KON .EQ. 0) THEN 342 343 IF (SPONLY .EQ. 0) THEN 344 CALL BLDNEWHS_SP 345 ELSE 346 CALL BLDNEWHS 347 ENDIF 348 349 ELSE 350 351 CALL KBLDNEWH 352 353 ENDIF 354 355 ! 356 ! If we're starting from a restart file, we need to modify H such 357 ! that it agrees with the density matrix elements read from file 358 ! 359 360 IF (RESTART .EQ. 1) CALL IFRESTART 361 362 ! 363 ! See whether we need spin-dependence too 364 ! 365 366 IF (SPINON .EQ. 1) THEN 367 CALL GETDELTASPIN 368 CALL BLDSPINH 369 ENDIF 370 371 IF (CONTROL .EQ. 1) THEN 372 CALL ALLOCATEDIAG 373 ELSEIF (CONTROL .EQ. 2 .OR. CONTROL .EQ. 4 .OR. CONTROL .EQ. 5) THEN 374 CALL ALLOCATEPURE 375 ELSEIF (CONTROL .EQ. 3) THEN 376 CALL FERMIALLOCATE 377 ENDIF 378 379 ! ELSE 380 ! CALL ERRORS("latte_lib","Attemting to perform multiple single point calculations. & 381 ! & MDON .EQ. 0 .AND. RELAXME .EQ. 0 can be done & 382 ! &for only one geometry (A single call to the library). Please reduce & 383 ! &the number of steps (md of relaxation) at the host code") 384 ! EXISTERROR_INOUT = EXISTERROR 385 ! RETURN 386 ! ENDIF 387 388 IF (CONTROL .EQ. 5) THEN 389 390 CALL GERSHGORIN 391 CALL SP2FERMIINIT 392 393 ENDIF 394 395 IF (ELECTRO .EQ. 0) CALL QNEUTRAL(0,1) ! Local charge neutrality 396 397 IF (ELECTRO .EQ. 1) CALL QCONSISTENCY(0,1) ! Self-consistent charges 398 399 ! We have to build our NKTOT complex H matrices and compute the 400 ! self consistent density matrix 401 402 ! Tr[rho dH/dR], Pulay force, and Tr[rho H] need to de-orthogonalized rho 403 404 IF (KON .EQ. 1) CALL KGETDOS 405 406 ! OPEN(UNIT=31, STATUS="UNKNOWN", FILE="myrho.dat") 407 408 ! DO I = 1, HDIM 409 ! DO J = 1,HDIM 410 411 ! IF (ABS(BO(J,I)) .GT. 1.0D-5) WRITE(31,99) I, J 412 413 ! ENDDO 414 ! ENDDO 415 416 !99 FORMAT(2I9) 417 !CLOSE(31) 418 419 IF (DEBUGON .EQ. 1 .AND. SPINON .EQ. 0 .AND. KON .EQ. 0) THEN 420 421 PRINT*, "Caution - you're writing to file the density matrix!" 422 423 OPEN(UNIT=31, STATUS="UNKNOWN", FILE="myrho.dat") 424 425 DO I = 1, HDIM 426 WRITE(31,10) (BO(I,J), J = 1, HDIM) 427 ENDDO 428 429 CLOSE(31) 430 43110 FORMAT(100G18.8) 432 433 ENDIF 434 435 FTOT = ZERO 436 437 IF (COMPFORCE .EQ. 1) THEN 438 439 IF (KON .EQ. 0) THEN 440 441 IF (SPONLY .EQ. 0) THEN 442 CALL GRADHSP 443 ELSE 444 CALL GRADH 445 ENDIF 446 447 ELSE 448 CALL KGRADH 449 ENDIF 450 451 FTOT = TWO * F 452 453 ENDIF 454 455 EREP = ZERO 456 IF (PPOTON .EQ. 1) THEN 457 CALL PAIRPOT 458 FTOT = FTOT + FPP 459 ENDIF 460 461 IF (PPOTON .EQ. 2) THEN 462 CALL PAIRPOTTAB 463 FTOT = FTOT + FPP 464 ENDIF 465 466 IF (PPOTON .EQ. 3) THEN 467 CALL PAIRPOTSPLINE 468 FTOT = FTOT + FPP 469 ENDIF 470 471 IF (ELECTRO .EQ. 1) FTOT = FTOT + FCOUL 472 473 IF (BASISTYPE .EQ. "NONORTHO") THEN 474 475 IF (SPONLY .EQ. 0) THEN 476 ! s/sp orbitals only so we can use the analytic code 477 CALL FCOULNONO_SP 478 CALL PULAY_SP 479 IF (SPINON .EQ. 1) CALL FSPINNONO_SP 480 481 FTOT = FTOT - TWO*FPUL 482 483 IF (ELECTRO .EQ. 0) FTOT = FTOT + FSLCN 484 IF (ELECTRO .EQ. 1) FTOT = FTOT + FSCOUL 485 IF (SPINON .EQ. 1) FTOT = FTOT + FSSPIN 486 487 ELSE 488 ! Otherwise use the complex but general expansions Josh 489 ! Coe implemented 490 491 CALL PULAY 492 493 FTOT = FTOT + FPUL 494 495! CALL FCOULNONO 496! CALL PULAY 497! IF (SPINON .EQ. 1) CALL FSPINNONO 498 ENDIF 499 500! FTOT = FTOT - TWO*FPUL 501 502! IF (ELECTRO .EQ. 0) FTOT = FTOT + FSLCN 503! IF (ELECTRO .EQ. 1) FTOT = FTOT + FSCOUL 504! IF (SPINON .EQ. 1) FTOT = FTOT + FSSPIN 505 506 507 !FTOT = FTOT - TWO*FPUL + FSCOUL 508 509 !IF (SPINON .EQ. 1) FTOT = FTOT + FSSPIN 510 511 ENDIF 512 513 CALL TOTENG 514 515 ECOUL = ZERO 516 IF (ELECTRO .EQ. 1) CALL GETCOULE 517 518 ESPIN = ZERO 519 IF (SPINON .EQ. 1) CALL GETSPINE 520 521 IF (CONTROL .NE. 1 .AND. CONTROL .NE. 2 .AND. KBT .GT. 0.000001 ) THEN 522 523 ! We get the entropy automatically when using diagonalization. 524 ! This is only required when employing the recursive expansion 525 ! of the Fermi-operator at finite electronic temperature 526 527 CALL ENTROPY 528 529 ENDIF 530 531 CALL WRTRESTART(0) 532 533 IF (CONTROL .EQ. 1) THEN 534 ! CALL DEALLOCATEDIAG 535 ELSEIF (CONTROL .EQ. 2 .OR. CONTROL .EQ. 4 .OR. CONTROL .EQ. 5) THEN 536 CALL DEALLOCATEPURE 537 ELSEIF (CONTROL .EQ. 3) THEN 538 CALL FERMIDEALLOCATE 539 ENDIF 540 541 ! 542 ! Stop the clocks 543 ! 544 545 TX = STOP_TIMER(LATTE_TIMER) 546 547 548#ifndef FCIDxlf 549 CALL DTIME(TARRAY, RESULT) 550#endif 551 552 CALL SYSTEM_CLOCK(STOP_CLOCK, CLOCK_RATE, CLOCK_MAX) 553 554 CALL GETPRESSURE 555 556 ! WRITE(*,*) "Force", FPP(1,1), FPP(2,1), FPP(3,1) 557 ! PRINT*, "PCHECK ", (1.0/3.0)*(VIRBOND(1)+VIRBOND(2) + VIRBOND(3)), & 558 ! (1.0/3.0)*(VIRCOUL(1)+VIRCOUL(2) + VIRCOUL(3)), & 559 ! (1.0/3.0)*(VIRPAIR(1)+VIRPAIR(2) + VIRPAIR(3)), & 560 ! (1.0/3.0)*(VIRPUL(1)+VIRPUL(2) + VIRPUL(3)), & 561 ! (1.0/3.0)*(VIRSCOUL(1)+VIRSCOUL(2) + VIRSCOUL(3)) 562 563#ifdef DBCSR_ON 564 565 IF (CONTROL .EQ. 2 .AND. SPARSEON .EQ. 1 .AND. MYNODE .EQ. 0) THEN 566 567#endif 568 569 IF (MYID .EQ. 0) THEN 570 CALL FITTINGOUTPUT(0) 571 CALL SUMMARY 572 573 ! IF (SPINON .EQ. 0) CALL NORMS 574 575 PRINT*, "# System time = ", TARRAY(1) 576 PRINT*, "# Wall time = ", FLOAT(STOP_CLOCK - START_CLOCK)/FLOAT(CLOCK_RATE) 577 PRINT*, "# Wall time per SCF =", & 578 FLOAT(STOP_CLOCK - START_CLOCK)/(FLOAT(CLOCK_RATE)*FLOAT(NUMSCF)) 579 ! PRINT*, HDIM, FLOAT(STOP_CLOCK - START_CLOCK)/FLOAT(CLOCK_RATE) 580 TX = TIMER_RESULTS() 581 PRINT*, "# NUMSCF = ", NUMSCF 582 583 ENDIF 584#ifdef DBCSR_ON 585 586 ENDIF 587 588#endif 589 590 ! CALL ASSESSOCC 591 592 IF (ELECTRO .EQ. 1) CALL DEALLOCATECOULOMB 593 594 IF (BASISTYPE .EQ. "NONORTHO") CALL DEALLOCATENONO 595 596 CALL DEALLOCATENEBARRAYS 597 598 CALL DEALLOCATEALL 599 600 LIBINIT = .FALSE. 601 602 RETURN 603 604 ELSEIF (MDON .EQ. 1 .AND. RELAXME .EQ. 0 .AND. MAXITER_IN < 0 ) THEN 605 606 IF(VERBOSE >= 1)WRITE(*,*)"Insie MDON= 1 and RELAXME= 0 ..." 607 608 DT = DT_IN ! Get the integration step from the hosting code. 609 610 V = VEL_IN/1000.0d0 !Convert from Ang/ps to Ang/fs 611 612 !Control for implicit geometry optimization. 613 !This will need to be replaced by a proper flag. 614 IF (DT_IN == 0) THEN 615 IF (VERBOSE >= 1) WRITE(*,*)"NOTE: DT = 0 => FULLQCONV = 1" 616 IF (VERBOSE >= 1) WRITE(*,*)"NOTE: DT = 0 => MDMIX = QMIX" 617 FULLQCONV = 1 618 MDMIX = QMIX 619 FLUSH(6) 620 ENDIF 621 622 IF (LIBCALLS == 0) THEN 623 624 IF (VERBOSE >= 1)WRITE(*,*)"Allocating nonorthogonal arrays ..." 625 IF (BASISTYPE .EQ. "NONORTHO") CALL ALLOCATENONO 626 627 IF (VERBOSE >= 1)WRITE(*,*)"Allocating XLBO arrays ..." 628 IF (XBOON .EQ. 1) CALL ALLOCATEXBO 629 630 IF (VERBOSE >= 1)WRITE(*,*)"Allocating COULOMB arrays ..." 631 IF (ELECTRO .EQ. 1) THEN 632 CALL ALLOCATECOULOMB 633 CALL INITCOULOMB 634 ENDIF 635 636 ! Start the timers 637 IF (VERBOSE >= 1)WRITE(*,*)"Starting timers ..." 638 CALL SYSTEM_CLOCK(START_CLOCK, CLOCK_RATE, CLOCK_MAX) 639 640#ifndef FCIDxlf 641 CALL DTIME(TARRAY, RESULT) 642#endif 643 644 IF (VERBOSE >= 1)WRITE(*,*)"Setting up TBMD ..." 645 CALL SETUPTBMD(NEWSYSTEM) 646 647 FLUSH(6) 648 649 ELSEIF (LIBCALLS > 0 .AND. RESTARTLIB == 0) THEN 650 651 DBOX = SQRT((BOX(1,1)-BOX_OLD(1,1))**2 + & 652 & (BOX(2,2)-BOX_OLD(2,2))**2 + (BOX(3,3)-BOX_OLD(3,3))**2) 653 654 ! Reinitializing Coulombic contribution if Kpoints are used 655 IF ((ELECTRO .EQ. 1 .AND. ELECMETH .EQ. 0) .OR. DBOX .GT. 0.0d0) THEN 656 CALL INITCOULOMB 657 ENDIF 658 659 IF (MOD(LIBCALLS, UDNEIGH) .EQ. 0) THEN 660 !If box is changing 661 IF (DBOX .GT. 0.0d0) THEN 662 CALL NEBLISTS(0) 663 CALL INITCOULOMB !If the box is changing we need to recompute the kspace list 664 ELSE 665 CALL NEBLISTS(1) 666 ENDIF 667 ENDIF 668 669 BOX_OLD = BOX 670 671 ENDIF 672 673 IF (QITER .NE. 0) THEN 674 ECOUL = ZERO 675 IF (ELECTRO .EQ. 1) CALL GETCOULE 676 ENDIF 677 678 IF(VERBOSE >= 1) WRITE(*,*)"LIBCALLS",LIBCALLS 679 680 IF(VERBOSE >= 1) MLSI = TIME_MLS() 681 IF(LIBCALLS > 0) CALL GETMDF(1, LIBCALLS) 682 IF(VERBOSE >= 1) WRITE(*,*)"Time for GETMDF =", TIME_MLS()-MLSI 683 684 CALL TOTENG 685 686 ! For the 0 SCF MD the coulomb energy is calculated in GETMDF 687 688 IF (PPOTON .EQ. 1) THEN 689 CALL PAIRPOT 690 ELSEIF (PPOTON .EQ. 2) THEN 691 CALL PAIRPOTTAB 692 ELSEIF (PPOTON .EQ. 3) THEN 693 CALL PAIRPOTSPLINE 694 ENDIF 695 696 IF (QITER .NE. 0) THEN 697 ECOUL = ZERO 698 IF (ELECTRO .EQ. 1) CALL GETCOULE 699 ENDIF 700 701 ESPIN = ZERO 702 IF (SPINON .EQ. 1) CALL GETSPINE 703 704 IF (CONTROL .NE. 1 .AND. CONTROL .NE. 2 .AND. KBT .GT. 0.000001 ) THEN 705 706 ! Only required when using the recursive expansion of the Fermi operator 707 708 ! 2/26/13 709 ! The entropy is now calculated when we get the density 710 ! matrix in the spin polarized case with diagonalization, 711 ! as it should be... 712 713 CALL ENTROPY 714 715 ENDIF 716 717 IF(VERBOSE >= 1) WRITE(*,*)"Energy Components (TRRHOH, EREP, ENTE, ECOUL)",TRRHOH, EREP, ENTE, ECOUL 718 719 IF (FREEZE .EQ. 1) CALL FREEZE_ATOMS(FTOT,V) 720 721 IF(MAXVAL(FTOT_OUT) .NE. 0.0d0)THEN 722 IF(VERBOSE >= 1) WRITE(*,*)"Adding force components and energies from applicacion code ..." 723 IF(VERBOSE >= 1) WRITE(*,*)"APPCODE,LATTE",VENERG,TRRHOH + EREP - ENTE - ECOUL + ESPIN 724 VENERG = TRRHOH + EREP - ENTE - ECOUL + ESPIN 725 FTOT_OUT = FTOT_OUT + FTOT 726 ELSE 727 VENERG = TRRHOH + EREP - ENTE - ECOUL + ESPIN 728 FTOT_OUT = FTOT 729 ENDIF 730 731 732 ! Get the seccond virial coefficient to pass it to the application program 733 IF (ELECTRO .EQ. 0) VIRCOUL = ZERO 734 VIRIAL = VIRBOND + VIRPAIR + VIRCOUL 735 736 IF (SPINON .EQ. 1) VIRIAL = VIRIAL + VIRSSPIN 737 738 IF (BASISTYPE .EQ. "NONORTHO") THEN 739 VIRIAL = VIRIAL - VIRPUL + VIRSCOUL 740 ENDIF 741 742 ! CALL GETPRESSURE 743 744 VIRIAL_INOUT = -VIRIAL 745 746 LIBINIT = .TRUE. 747 NEWSYSTEM = 0 !Setting newsystem back to 0. 748 749#ifdef PROGRESSON 750 IF(MOD(LIBCALLS,WRTFREQ) == 0)THEN 751 IF(VERBOSE >= 1) THEN 752 WRITE(*,*)"Writing trajectory into trajectory.pdb ..." 753 SY%NATS = NATS 754 IF(.NOT. ALLOCATED(SY%COORDINATE))ALLOCATE(SY%COORDINATE(3,NATS)) 755 SY%COORDINATE = CR 756 SY%SYMBOL = ATELE 757 SY%LATTICE_VECTOR = BOX 758 SY%NET_CHARGE = DELTAQ 759 CALL PRG_WRITE_TRAJECTORY(SY,LIBCALLS,WRTFREQ,DT_IN,"trajectory","pdb") 760 761 WRITE(*,*)"Writing trajectory into trajectory.xyz ..." 762 IF(LIBCALLS .EQ. 0)THEN 763 OPEN(UNIT=20,FILE="trajectory.xyz",STATUS='unknown') 764 ELSE 765 OPEN(UNIT=20,FILE="trajectory.xyz",POSITION='append',STATUS='old') 766 ENDIF 767 !Extended xyz file. 768 WRITE(20,*)NATS 769 WRITE(20,*) 'Lattice="',BOX(1,1),BOX(1,2),BOX(1,3),& 770 &BOX(2,1),BOX(2,2),BOX(2,3),BOX(3,1),BOX(3,2),BOX(3,3),'"',& 771 &"Properties=species:S:1:pos:R:3:vel:R:3:for:R:3:cha:R:1 Time=",LIBCALLS*DT_IN 772 DO I=1,NATS 773 WRITE(20,*)ATELE(I),CR(1,I),CR(2,I),CR(3,I),V(1,I),V(2,I),V(3,I),& 774 &FTOT(1,I),FTOT(2,I),FTOT(3,I),-DELTAQ(I) 775 ENDDO 776 CLOSE(20) 777 ENDIF 778 ENDIF 779#else 780 IF(MOD(LIBCALLS,WRTFREQ) == 0)THEN 781 IF(VERBOSE >= 1) THEN 782 WRITE(*,*)"Writing trajectory into trajectory.xyz ..." 783 IF(LIBCALLS .EQ. 0)THEN 784 OPEN(UNIT=20,FILE="trajectory.xyz",STATUS='unknown') 785 ELSE 786 OPEN(UNIT=20,FILE="trajectory.xyz",POSITION='append',STATUS='old') 787 ENDIF 788 !Extended xyz file. 789 WRITE(20,*)NATS 790 WRITE(20,*) 'Lattice="',BOX(1,1),BOX(1,2),BOX(1,3),& 791 &BOX(2,1),BOX(2,2),BOX(2,3),BOX(3,1),BOX(3,2),BOX(3,3),'"',& 792 &"Properties=species:S:1:pos:R:3:vel:R:3:for:R:3:cha:R:1 Time=",LIBCALLS*DT_IN 793 DO I=1,NATS 794 WRITE(20,*)ATELE(I),CR(1,I),CR(2,I),CR(3,I),V(1,I),V(2,I),V(3,I),& 795 &FTOT(1,I),FTOT(2,I),FTOT(3,I),-DELTAQ(I) 796 ENDDO 797 CLOSE(20) 798 ENDIF 799 ENDIF 800#endif 801 802 IF(VERBOSE >= 1 .AND. CONTROL == 1 .AND. KON == 0)THEN 803 IF(SPINON == 0) THEN 804 HOMO = EVALS(FLOOR(BNDFIL*FLOAT(HDIM))) 805 LUMO = EVALS(FLOOR(BNDFIL*FLOAT(HDIM))+1) 806 WRITE(*,*)"HOMO=",HOMO, "LUMO=",LUMO 807 WRITE(*,*)"EGAP=",LUMO - HOMO 808 ELSE 809 HOMO = MAX(DOWNEVALS(FLOOR(BNDFIL*FLOAT(HDIM))),UPEVALS(FLOOR(BNDFIL*FLOAT(HDIM)))) 810 LUMO = MIN(DOWNEVALS(FLOOR(BNDFIL*FLOAT(HDIM))+1),UPEVALS(FLOOR(BNDFIL*FLOAT(HDIM))+1)) 811 WRITE(*,*)"HOMO=",HOMO, "LUMO=",LUMO 812 WRITE(*,*)"EGAP=",LUMO - HOMO 813 ENDIF 814 ENDIF 815 816 IF (MOD(LIBCALLS, RSFREQ) .EQ. 0)THEN 817 IF(VERBOSE >= 0) CALL WRTRESTARTLIB(LIBCALLS) 818 ENDIF 819 820 IF(RESTARTLIB == 1 .AND. LIBCALLS == 0)THEN 821 CALL READRESTARTLIB(LIBCALLS) 822 ENDIF 823 824#ifdef PROGRESSON 825 ! IF(VERBOSE >= 1)THEN 826 ! CALL GETDIPOLE(DIPOLEMAG,DIPOLEVECOUT=DIPOLEVECOUT) 827 ! WRITE(*,*)"Dipole Magnitude = DIPOLEMAG" 828 ! WRITE(*,*)"Dipole Vector:" 829 ! WRITE(*,*)DIPOLEVECOUT(1),DIPOLEVECOUT(2),DIPOLEVECOUT(3) 830 ! ENDIF 831#endif 832 833 FLUSH(6) !To force writing to file at every call 834 835 EXISTERROR_INOUT = EXISTERROR 836 837 RETURN 838 839 ELSEIF (MDON .EQ. 1 .AND. RELAXME .EQ. 0 .AND. MAXITER_IN >= 0) THEN 840 841 IF (BASISTYPE .EQ. "NONORTHO") CALL ALLOCATENONO 842 843 IF (XBOON .EQ. 1) CALL ALLOCATEXBO 844 845 IF (ELECTRO .EQ. 1) THEN 846 CALL ALLOCATECOULOMB 847 CALL INITCOULOMB 848 ENDIF 849 850 ! Start the timers 851 852 CALL SYSTEM_CLOCK(START_CLOCK, CLOCK_RATE, CLOCK_MAX) 853 854#ifndef FCIDxlf 855 CALL DTIME(TARRAY, RESULT) 856#endif 857 858 ! 859 ! Call TBMD 860 ! 861 862 CALL TBMD 863 864#ifdef MPI_ON 865 IF (PARREP .EQ. 1) CALL MPI_BARRIER (MPI_COMM_WORLD, IERR ) 866#endif 867 868 ! Stop the timers 869 870#ifndef FCIDxlf 871 CALL DTIME(TARRAY, RESULT) 872#endif 873 874 CALL SYSTEM_CLOCK(STOP_CLOCK, CLOCK_RATE, CLOCK_MAX) 875 876 CALL SUMMARY 877 878 IF (PBCON .EQ. 0) CLOSE(23) 879 880 IF (BASISTYPE .EQ. "NONORTHO") CALL DEALLOCATENONO 881 882 IF (XBOON .EQ. 1) CALL DEALLOCATEXBO 883 884 IF (ELECTRO .EQ. 1) CALL DEALLOCATECOULOMB 885 886 ! SYSTPURE = TARRAY(1) 887 ! WRITE(6,'("# System time for MD run = ", F12.2, " s")') SYSTPURE 888 WRITE(6,'("# Wall time for MD run = ", F12.2, " s")') & 889 FLOAT(STOP_CLOCK - START_CLOCK)/FLOAT(CLOCK_RATE) 890 891 892 ELSEIF (MDON .EQ. 0 .AND. RELAXME .EQ. 1) THEN 893 894 CALL ERRORS("latte_lib","This option was not tested for the & 895 &library version of LATTE: RELAXME= 1") 896 RETURN 897 898 CALL MSRELAX 899 900 ELSEIF (MDON .EQ. 0 .AND. RELAXME .EQ. 0 .AND. DOSFITON .EQ. 1) THEN 901 902 CALL ERRORS("latte_lib","This option was not tested for the & 903 &library version of LATTE: DOSFITON= 1") 904 RETURN 905 906 CALL SYSTEM_CLOCK(START_CLOCK, CLOCK_RATE, CLOCK_MAX) 907 908 CALL DOSFIT 909 910 CALL SYSTEM_CLOCK(STOP_CLOCK, CLOCK_RATE, CLOCK_MAX) 911 912 WRITE(6,'("# Wall time = ", F12.2, " s")') & 913 FLOAT(STOP_CLOCK - START_CLOCK)/FLOAT(CLOCK_RATE) 914 915 ELSEIF (MDON .EQ. 0 .AND. RELAXME .EQ. 0 .AND. DOSFITON .EQ. 2) THEN 916 917 CALL ERRORS("latte_lib","This option was not tested for the & 918 &library version of LATTE: DOSFITON= 3") 919 RETURN 920 921 CALL MOFIT 922 923 ELSEIF (MDON .EQ. 0 .AND. RELAXME .EQ. 0 .AND. DOSFITON .EQ. 3) THEN 924 925 CALL ERRORS("latte_lib","This option was not tested for the & 926 &library version of LATTE: DOSFITON= 3") 927 RETURN 928 929 CALL MOFITPLATO 930 931 ELSEIF (MDON .EQ. 0 .AND. RELAXME .EQ. 0 .AND. PPFITON .EQ. 1) THEN 932 933 CALL ERRORS("latte_lib","This option was not tested for the & 934 &library version of LATTE: PPFITON= 1") 935 RETURN 936 937 CALL PPFIT 938 939 ELSEIF (MDON .EQ. 0 .AND. RELAXME .EQ. 0 .AND. ALLFITON .EQ. 1) THEN 940 941 CALL ERRORS("latte_lib","This option was not tested for the & 942 &library version of LATTE: ALLFITON= 1") 943 RETURN 944 945 CALL ALLFIT 946 947 ELSE 948 949 CALL ERRORS("latte_lib","You can't have RELAXME = 1 and MDON = 1") 950 951 ENDIF 952 953#ifdef GPUON 954 955 CALL SHUTDOWN() 956 957#endif 958 959 960 CALL DEALLOCATEALL 961 962#ifdef DBCSR_ON 963 964 !ends mpi 965 966 IF (CONTROL .EQ. 2 .AND. SPARSEON .EQ. 1) CALL SHUTDOWN_DBCSR 967 968#endif 969 970 ! Done with timers 971 TX = SHUTDOWN_TIMER() 972 973#ifdef MPI_ON 974 CALL MPI_FINALIZE( IERR ) 975#endif 976 977 LIBINIT = .TRUE. 978 NEWSYSTEM = 0 !Setting newsystem back to 0. 979 EXISTERROR_INOUT = EXISTERROR 980 981 END SUBROUTINE LATTE 982 983END MODULE LATTE_LIB 984