aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorRahiel Kasim <rahielkasim@gmail.com>2016-02-18 00:26:07 +0100
committerRahiel Kasim <rahielkasim@gmail.com>2016-02-18 00:26:07 +0100
commitca07195c57a3a0923d03c70af44b846eb15a9bd3 (patch)
tree33e98e69a98c584fc65cb814b835bbdc44a947d3
parent8d581178ea4014ffc0a108eb7677ee06b89f0bae (diff)
calculate energy drift from equilibrium due to numeric integration
-rwxr-xr-xSource/md.f20
1 files changed, 19 insertions, 1 deletions
diff --git a/Source/md.f b/Source/md.f
index a5fa60d..781e629 100755
--- a/Source/md.f
+++ b/Source/md.f
@@ -24,7 +24,7 @@ c
LOGICAL scale, multi_on
DOUBLE PRECISION en, ent, vir, virt, enk, time, enpot, delt, tmax,
& enkt, tequil, temprsq
- DOUBLE PRECISION enpot2, vir2
+ DOUBLE PRECISION enpot2, en2, vir2, enk2
c --common blocks declaration:
INCLUDE 'parameter.inc'
INCLUDE 'conf.inc'
@@ -37,7 +37,11 @@ c --common blocks declaration:
INTEGER nlist2(npmax), list2(npmax, npmax)
DOUBLE PRECISION fx(NPMax), fy(NPMax), fz(NPMax)
DOUBLE PRECISION fx2(NPMax), fy2(NPMax), fz2(NPMax)
+ INTEGER E0_step
+ DOUBLE PRECISION E0, E_drift
+ E0_step = -1 ! steps since we recorded E_0 after equilibration
+ E_drift = 0.D0
WRITE (6, *) '**************** MC_NPT ***************'
c ---initialize system
@@ -77,6 +81,17 @@ c Verlet-list
END IF
en = enpot + enk
step = step + 1
+ IF (E0_step.eq.-1 .and. time.gt.tequil) THEN
+c --- Set initial energy after equilibration
+ CALL TOTERG(en2, vir2, enk2)
+ E0 = en2
+ E0_step = 0
+ ELSE IF (E0_step.gt.-1) THEN
+c --- Calculate energy drift
+ CALL TOTERG(en2, vir2, enk2)
+ E_drift = E_drift + abs((en2 - E0) / E0)
+ E0_step = E0_step + 1
+ END IF
IF (time.LT.tequil) THEN
IF (scale) THEN
IF (MOD(step,20).EQ.0) CALL VELOCS(temprsq)
@@ -99,6 +114,9 @@ c ---write intermediate configuration to file
IF (SAMP2) CALL SAMPLE2(2, delt)
WRITE (6, 99002) ent, virt
CALL STORE(21)
+c --- Report average energy drift from the initial energy
+ E_drift = E_drift / (E0_step - 1)
+ write (6, *) "Average Energy Drift: ", E_drift
STOP
99001 FORMAT (' Total pot. energy in. conf. : ', f12.5, /,