EVB_Project/evb_input at master · glass-w/EVB_Project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
*AADH 2 state EVB
*
!
!            Amine                             Amine
!              |                                 |
!              NH                                NH
!              |
!              H
!                                             H
!                                            /
!           O(D1)  O(D2)                     O(D1)  O(D2)
!           \     //                          \     //
!            \   //                            \   //
!              C3                               C3
!              |                                |
!             Asp                              Asp
!      state 1 (aadh_reactant.psf)      state 2 (aadh_prod_od1.psf) !-- these are the psf files

! nensem : number of EVB states
ensemble nensem 2

bomlev 0
prnlev 5
wrnlev 2

if ?ensemble .ne. 1 then
 echo "Test NOT performed."
 stop
endif

set node ?whoiam
set setupdir /Users/williamglass/Documents/TMCS/Bristol/Project/AADH_EVB
set setupdirOD1 /Users/williamglass/Documents/TMCS/Bristol/Project/AADH_EVB/OD1
set outputdir /Users/williamglass/Documents/TMCS/Bristol/Project/AADH_EVB/OD1_outputs

if @node .EQ. 0 set top state1
if @node .EQ. 1 set top state2

!
! READ IN TOPOLOGIES
!
ensemble open read unit 11 card name @setupdir/top_ttw@node.inp
read rtf card unit 11
ensemble close unit 11

!
! READ IN PARAMETERS
!

ensemble open read form unit 1 name @setupdir/par_ttw.inp
read parameter card unit 1
ensemble close unit 1

!
! READ IN PSF & COOR
!
if @node .eq. 0 then
read psf card name @setupdirOD1/aadh_reactant.psf
endif
if @node .eq. 1 then
read psf card name @setupdir/aadh_prod_od1.psf
endif

! You only need on coordinate file.  Eveything should have the same coordinate.
read coor card name @setupdirOD1/path1_b3lyp_geom1.crd
!
! CALCULATE DIPOLE
!
!  This is needed to allow the calculation of state averaged dipole moments.
!  Useful in debugging/analysis. As there is only one molecule then we can use
!  'sele all end'
!
! coor dipo sele all end

!
! CALCULATE NON-EVB ENERGY
!
!  This is here for comparison with EVB energy calculation
!
energy

!
! SET UP EVB CALCULATION
!
!  Coupling parameters for 2-D Gaussian + a constant coupling
!
! coupling of zero will give a EVB matrix with the energies of the states as the
! eigenvalues.

set CONS 0.0

!
! INITIALISE EVB CALCULATION
!
!  EVB hamiltonian matrix, H = :
!
!  state 1: |[0,0] [0,1]|
!  state 2: |[1,0] [1,1]|
!
!  Only modify upper triangular part as lower triangular part is filled in
!  automatically.
!
!   V0 and V1 and the potential energies of state 1 and 2 respectively.
!
!   EVBS : save unit for dipole output
!   SHFT 0 <VALUE>                      : H[0,0] <- V0 + <VALUE>
!   COUP 0 1 CONS <VALUE>               : H[0,1] <- H[0,1] + <VALUE>
!   COUP 0 1 GAUS TWOD <2D GAUSSIAN>    : H[0,1] <- H[0,1] + <2D GAUSSIAN>

ensemble evb
       COUP 0 1 CONS @CONS
envbi

energy

!
! RUN RELAXED SCAN
!
!ensemble open write card unit 21 name @outputdir/reaction.energy.@node

!set atom1 mal 1 O1
!set atom2 mal 1 H
!set atom3 mal 1 O2

!set inc 0.1
!set v -1.0
!set n 1

!label loop
!       skip none

!       resd reset kval 2000.0 rval @v -
!       1.0 @atom1 @atom2 -1.0 @atom2 @atom3

!       mini abnr nstep 200 nprint 100

!       skip resd
!       energy
!       write title unit 21
!       * @v ?ener
!       *

!incr v by @inc
!incr n by 1
!if n .le. 21 goto loop

!ensemble close unit 21

stop