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 22SUBROUTINE NVTRESCALE 23 24 USE CONSTANTS_MOD 25 USE SETUPARRAY 26 USE MDARRAY 27 USE MYPRECISION 28 29 IMPLICIT NONE 30 31 INTEGER :: I, J, K 32 REAL(LATTEPREC) :: PREFACTOR, PREPREFACT 33 REAL(LATTEPREC) :: CHI, TMPTEMP, MYTEMP, MOMENTUM(3) 34 REAL(LATTEPREC) :: MAXCHI = 1.05D0 35 IF (EXISTERROR) RETURN 36 37 PREPREFACT = HALF*F2V*DT 38 39 TMPTEMP = ZERO 40 41 DO I = 1, AVEPER 42 TMPTEMP = TMPTEMP + THIST(I) 43 ENDDO 44 45 MYTEMP = TMPTEMP/FLOAT(AVEPER) 46 47 ! CHI = SQRT(TTARGET/MYTEMP) 48 49 ! CHI = MIN(CHI, MAXCHI) 50 CHI = MIN(SQRT(TTARGET/MYTEMP), MAXCHI) 51 52 ! PRINT*, "MYTEMP = ", MYTEMP, "CHI = ", CHI 53 54 DO I = 1, NATS 55 56 PREFACTOR = PREPREFACT/MASS(ELEMPOINTER(I)) 57 58 ! 59 ! Half timestep advance in V with velocity rescale 60 ! 61 62 V(1,I) = CHI*V(1,I) + PREFACTOR*FTOT(1,I) 63 V(2,I) = CHI*V(2,I) + PREFACTOR*FTOT(2,I) 64 V(3,I) = CHI*V(3,I) + PREFACTOR*FTOT(3,I) 65 66 ENDDO 67 68 ! 69 ! Whole timestep advance in positions 70 ! 71 72 CR = CR + DT*V 73 74 ! 75 ! Get new force to complete advance in V 76 ! 77 78 CALL GETMDF(1, 100) 79 80 ! 81 ! Now finish advancing V with F(t + dt) 82 ! 83 84 DO I = 1, NATS 85 86 PREFACTOR = PREPREFACT/MASS(ELEMPOINTER(I)) 87 88 V(1,I) = V(1,I) + PREFACTOR*FTOT(1,I) 89 V(2,I) = V(2,I) + PREFACTOR*FTOT(2,I) 90 V(3,I) = V(3,I) + PREFACTOR*FTOT(3,I) 91 92 ENDDO 93 94 ! Remove drift 95 96 MOMENTUM = ZERO 97 98 DO I = 1, NATS 99 MOMENTUM(1) = MOMENTUM(1) + V(1,I)*MASS(ELEMPOINTER(I)) 100 MOMENTUM(2) = MOMENTUM(2) + V(2,I)*MASS(ELEMPOINTER(I)) 101 MOMENTUM(3) = MOMENTUM(3) + V(3,I)*MASS(ELEMPOINTER(I)) 102 ENDDO 103 104 MOMENTUM = MOMENTUM/SUMMASS 105 106 DO I = 1, NATS 107 V(1,I) = V(1,I) - MOMENTUM(1) 108 V(2,I) = V(2,I) - MOMENTUM(2) 109 V(3,I) = V(3,I) - MOMENTUM(3) 110 ENDDO 111 112 RETURN 113 114END SUBROUTINE NVTRESCALE 115 116