EHMACC, a QCPE program.
We have seamlessly incorporated the code into ours. Hence, unlike
the CNDO/ET code that separates the electronic structure part from
the electronic dynamics part, the EH/ET code can do both in a single
package.
The main subroutines are summarized here:
main.f: The main program.
readh(tau,natom,ndim,nstp,tmd,st,ht,ft): Reads the time
series of PDB files, calls smat and hmat,
sets up the Hamiltonian series ht and the overlap series
st, and stores the trajectory into the common block
traj.
getpos(natom,nstp): Gets the atomic positions from
the PDB file series.
getpar(natom,ndimen,nstp): Gets the atomic parameters
for EH calculation, from the block data
that defines the atomic parameters
(i.e., the number of valence electrons, the valence shell ionization
potentials, the Slater type orbital exponents, and the weighting
coefficients), which were taken from the table collected by S.
Alvarez, Dept. de Quimica Inorganica, Univ. de Barcelona, Spain
smat(natom2,natom,ndim): Constructs the overlap matrix.
It calls the following subroutines:
movlap(natom2,natom,ndim): calculates the overlaps;
mov(sigma,pi,delta,phi,i,j,rr,na,nb,la,lb,nn):
calculates the overlap component independent of the angle
between the atoms;
abfns(a,b): calculate the A, B functions;
lovlap(strad,a,b): called from mov;
hmat(natom,ndim): Constructs the Huckel
Hamiltonians from overlaps.
difov(natom,ndim,mtraj,ntraj): Computes (atomic)
advanced or retarded
overlap matrices. mtraj and ntraj stand for
two points along the trajectory.
friction(natom,ndim,nstp,tmd,ft): Computes the
nonadiabatic electronic
friction matrix (the K matrix) using third-order differentiation in
the AO basis.
fribo(ndim,nstp,tmd,ft): Transforms the K matrix
to the NBO basis.
part1(natom,ndim,nstp,tmd,st,ht,st2,ht2): Diagonalizes
the Hamiltonian time series step by step. The main entry to the static
calculation part.
denmat(natom,ndim,n,st,ht): Calculates the density matrix,
the molecular orbital populations, and the gross populations at atoms.
grn(istep,natom,ndim,e,green): Calculates the Green
functions, e is the energy argument for the Green function
export green.
caltda(natom,ndim,e,tda,hdd,haa): Calculates the
superexchange coupling TDA using Larsson's equation,
TDD and TAA are also output.
subdiag(natom,ndim,t,nstp,aring,nring): Constructs
orthogonal matrix T for transformation from AO's to natural hybrid bond
orbitals, using input density matrix. Maximum number of aromatic
rings allowed = 10.
diagm(ellow,elup,invs,all,eonly): Solves the secular
equation for the general Hermitian eigenvalue problem. For systems with
inversion (i.e. that reduce to a real symmetric matrix), real versions
of the routines are used. For all=.true.,
all eigenvalues and eigenvectors are determined.
For all=.false., only those between ellow
and elup are determined. If eonly=.true.,
then no eigenvectors are determined (only has an effect if
all=.false.). This is my favorite diagonalization subroutine
in FORTRAN. It never disappoints me. In the NBO part of the EH/ET code,
we don't use JACVEC any longer.
load(isel,iat1,iat2,dblk,sblk,nb): Zeroes the 8x8 matrices
dblk and sblk,
loads in atomic blocks of density matrix (when isel=1),
Hamiltonian matrix (when isel=2), and overlap matrix for the
atoms iat1 and iat2.
ulll(natom,ndim): Finds the orbit index upper and lower
bounds ul and ll.
deplet(borb,occ,nb,iat1,iat2): Depletes from the density
matrix the contributions from lone pairs.
STASH(BORB,OCC,IBD,IAT1,IAT2,IAT3): Seperates bond
orbital BORB into polarization coefficients
(stored in POL) and renormalized
hybrids (stored in Q),
keeps record of the number INO(NA) of hybrids
accumulated for each atom NA.
fourth(ibd,iring): Builds the fourth vector perpendicular
to the three known sp2 hydrids for all the carbon atoms
at the iring-th aromatic ring, and stores it in the
Q matrix.
reform(t,natom,ndim,nring): Builds the final orthogonal
transformation matrix T from polarization coefficients
(POL) and atomic hybrids (stored in Q,
identified in IATHY),
for bonds and their antibond counterparts. Note that the subblock
diagonalization method identifies only the bonding orbitals,
the antibonds are then
constructed from the bonds by use of the orthogonality conditions.
Therefore, in the NBO, all antibonds are strictly orthogonal to their
bond counterparts. This is different from the fact that the natural
bonds constructed by the subblock diagonaliztion method are only
marginally orthogonal to each other.
reformaro(t,ndim,nring): Establishes the
aromatic rings part of the T matrix.
anlyze(t,ndim): Finally analyzes the T matrix.
It calls:
htype(h,coef,pow,theta,phi): returns coefficients,
p-characters, and directions of hybrids (h(i),i=1,4)
from the transformation matrix. pow=0,1,2,...
for s, sp, sp2, ... hybrids
(pow=100 for pure p orbital),
coef = coefficient of hybrid,
theta,phi = polar and azimuthal angles of directed hybrid.
Normally pow is not strictly an integer. The bigger it is,
the more p character the hybrid orbital contains.
findsign(ndim,t,signname): Remembers the orientations
of all natural bonds at the initial time.
compsign(ndim,t): Fixes the signs for the natural hybrids
to avoid the random orientation flipping problem. A comparison phase is
obtained from a reference structure. A sign array which carries the
orientational information of the bonds is compared with the signs
stored in the comparison array. Bond orientations are corrected if
any disagreement. (read a detailed
explanation.)
compsignaro(ndim,t) is a separate subroutine special for
dealing with aromatic rings.
trans(ndim,t,n,st,ht,aring,nring): Transforms overlap,
Hamiltonian, and density matrices to the NBO basis.
subhdba(natom,ndim): Partitions Hamiltonian into
submatrices for donor, acceptor and bridge.
bspace(ndim,nstp): Diagonalizes the bridge Hamiltonian,
finds out the LUMO and HOMO of the bridge electronic structure.
faster(ndim,nstp,diff): Faster calculation of
propogation in static case, using the analytical formula of Green functions
in the time domain GDA(t).
This subroutine can be used to check the accuracy of the integration.
part2(tau,natom,ndim,nstp,tmd,st,ht,ft,st2,ht2):
The main entry for the electronic dynamics part.
hint(ndim,nstp,tmd,st,ht,ft,vtr,vti,cs,ch,cvr,cvi):
Prepares interpolation coefficients from the time series
for the overlap matrices, the Huckel matrices, together with
S-1(t)H(t) (vtr) and
S-1(t)K(t) (vti).
icsccu(x,y,nx,c,ic,ier): The cubic spline interpolation
subroutine. x,y are the input arrays, the output array
c stores the spline coefficients.
dynamics(tau,ndim,nstp,tmd,st,ht,vtr,vti,cs,ch,cvr,cvi):
Solves the time-dependent Schroedinger equation. Note that in nonorthogonal
basis, the site occupancies are not simply the module of the wave functions,
they should be a multiplication of the overlap matrix with the module.
dyntsa(tau,nstp,tmd,st2,ht2): Does the two-state dynamics.
strassen(ndim,nhaf,a,b,c): Matrix multiplication: The
Strassen method. The Strassen method is faster than gaxpy.
Since not every machine has BLAS library (hence dgemm),
it is an acceptable cross-platform solution for matrix multiplication.
predictor(tau,n,psi) and
corrector(tau,n,psi,time,nstp,tmd,vtr,vti,cvr,cvi):
Solves the second order differential equations using Gear's
predictor-corrector method.
rk4(tau,n,psi,time,nstp,tmd,vtr,vti,cvr,cvi):
Solves the second order differential equations using the classical
fourth order Runge-Kutta method.
HP-UX version and
Linux version
Remedy the gap problem
Usually the gap predicted by the EH method is about 4 eV (in contrast
to approximately 12 eV by the CNDO method), immediately putting itself
under criticism. We found, by chance, that the underestimation arises
the failure of treating the pi bonds. If one analyzes the states
near the LUMO and HOMO, one would find they actually are nothing but
the pi bonds. Here is an example:
1 BD 1.9859 -18.3210 C 1 H 2
2 BD 1.9963 -18.3619 C 1 H 3
3 BD 2.0075 -18.3691 C 1 H 4
.
.
.
111 LP 1.9962 -22.8611 O 6
112 LP 1.9687 -14.9379 O 6
113 LP 1.6838 -13.4116 N 7
.
.
.
139 BD* -.0113 8.6579 C 1 H 2
140 BD* -.0109 8.7869 C 1 H 3
141 BD* -.0111 8.8078 C 1 H 4
142 BD* .0036 14.3201 C 1 C 5
143 BD* .0087 7.6065 C 5 O 6
144 BD* .3747 -10.3492 C 5 O 6--> 1.6508
145 BD* .0067 10.1146 C 5 N 7
146 BD* -.0105 7.6999 N 7 H 8
147 BD* -.0163 6.5060 N 7 C 9
148 BD* .0014 11.2253 C 9 H 10
149 BD* -.0210 12.1555 C 9 C 11
150 BD* .0131 14.1056 C 9 C 22
151 BD* -.0054 8.7886 C 11 H 12
152 BD* -.0053 8.5700 C 11 H 13
153 BD* -.0162 10.7885 C 11 C 14
154 BD* .0038 9.0592 C 14 H 15
155 BD* -.0001 8.6885 C 14 H 16
156 BD* .0348 -8.2898 C 14 O 17--> 3.7102
157 BD* -.0041 -8.2409 O 17 C 18--> 3.7591
-----156 and 157 have p-characters of about 4
and 11, respectively. see the log file.
They are not C=O pi bonds. (We replaced S
with O because we couldn't treat S properly.)
158 BD* -.0015 9.0922 C 18 H 19
159 BD* -.0051 8.8397 C 18 H 20
160 BD* -.0025 8.8345 C 18 H 21
161 BD* .0077 11.0615 C 22 O 23
162 BD* .3804 -10.1526 C 22 O 23--> 1.8474
163 BD* -.0014 16.5099 C 22 N 24
164 BD* -.0072 7.6252 N 24 H 25
165 BD* -.0198 6.1701 N 24 C 26
166 BD* .0037 10.9273 C 26 H 27
167 BD* -.0105 9.7218 C 26 C 28
168 BD* .0143 11.1053 C 26 C 43
169 BD* -.0003 8.8298 C 28 H 29
170 BD* -.0037 8.6066 C 28 H 30
171 BD* -.0222 11.2918 C 28 C 31
172 BD* -.0010 9.0296 C 31 H 32
173 BD* -.0048 8.1757 C 31 H 33
174 BD* -.0263 16.8034 C 31 C 34
175 BD* -.0007 9.0944 C 34 H 35
176 BD* -.0042 7.9843 C 34 H 36
177 BD* -.0194 10.6734 C 34 C 37
178 BD* .0076 8.6021 C 37 H 38
179 BD* .0064 9.3603 C 37 H 39
180 BD* -.0186 5.1731 C 37 N 40
181 BD* -.0186 4.2648 N 40 H 41
182 BD* -.0128 5.0667 N 40 H 42
183 BD* .0083 9.0537 C 43 O 44
184 BD* .4143 -10.3166 C 43 O 44--> 1.6834
185 BD* -.0001 19.6593 C 43 N 45
186 BD* -.0159 7.4730 N 45 H 46
187 BD* -.0244 5.4937 N 45 C 47
188 BD* -.0023 11.2935 C 47 H 48
189 BD* -.0037 10.9374 C 47 H 49
190 BD* .0043 16.9445 C 47 C 50
191 BD* .0065 7.2201 C 50 O 51
192 BD* .3822 -10.3748 C 50 O 51--> 1.6252
193 BD* .0028 12.0287 C 50 N 52
194 BD* -.0127 7.5023 N 52 H 53
195 BD* -.0158 6.0018 N 52 C 54
196 BD* -.0012 10.6199 C 54 H 55
197 BD* .0001 10.7547 C 54 C 56
198 BD* .0052 17.1782 C 54 C 64
199 BD* .0134 9.4255 C 56 H 57
200 BD* -.0083 -.0686 C 56 O 58
201 BD* -.0094 7.4706 C 56 C 60
202 BD* -.0141 1.4505 O 58 H 59
203 BD* -.0074 9.0127 C 60 H 61
204 BD* -.0083 8.9262 C 60 H 62
205 BD* -.0089 8.5174 C 60 H 63
206 BD* .0066 9.4715 C 64 O 65
207 BD* .4263 -10.2937 C 64 O 65--> 1.7063
208 BD* .0017 18.1604 C 64 N 66
209 BD* -.0169 7.5421 N 66 H 67
210 BD* -.0145 5.4272 N 66 C 68
211 BD* .0001 10.7664 C 68 H 69
212 BD* -.0137 9.5187 C 68 C 70
213 BD* .0087 16.3618 C 68 C 83
214 BD* -.0029 8.8512 C 70 H 71
215 BD* -.0027 8.6398 C 70 H 72
216 BD* -.0182 11.4511 C 70 C 73
217 BD* -.0006 8.4893 C 73 H 74
218 BD* -.0191 5.8586 C 73 C 75
219 BD* -.0239 12.8016 C 73 C 79
220 BD* -.0093 8.4317 C 75 H 76
221 BD* -.0075 8.9759 C 75 H 77
222 BD* -.0030 9.5937 C 75 H 78
223 BD* -.0053 9.0430 C 79 H 80
224 BD* -.0064 8.7029 C 79 H 81
225 BD* -.0073 8.6264 C 79 H 82
226 BD* .0082 10.0190 C 83 O 84
227 BD* .3862 -10.2289 C 83 O 84--> 1.7711
228 BD* .0084 14.4269 C 83 N 85
229 BD* -.0158 7.5028 N 85 H 86
230 BD* -.0151 4.2080 N 85 C 87
231 BD* -.0014 10.8179 C 87 H 88
232 BD* -.0098 14.0006 C 87 C 89
233 BD* .0066 16.3324 C 87 C 97
234 BD* .0119 9.3137 C 89 H 90
235 BD* -.0022 -1.1148 C 89 O 91
236 BD* -.0078 8.8405 C 89 C 93
237 BD* -.0042 1.2832 O 91 H 92
238 BD* -.0079 8.7263 C 93 H 94
239 BD* -.0111 8.2585 C 93 H 95
240 BD* -.0066 8.8269 C 93 H 96
241 BD* .0094 11.5446 C 97 O 98
242 BD* .3371 -10.0943 C 97 O 98--> 1.9057
243 BD* .0101 11.6049 C 97 N 99
244 BD* -.0124 7.7769 N 99 H100
245 BD* -.0229 4.9143 N 99 C101
246 BD* .0010 9.1270 C101 H102
247 BD* -.0013 9.1258 C101 H103
248 BD* -.0042 9.1670 C101 H104
--------------------------------------------
minimum bond occupancy= 1.6555 at bond 130
maximum anti occupancy= .4263 at bond 207
------These two numbers mean that the electronic
wave functions of bonds are really very
localized. The Huckel approximation, as
the zero order approximation of the more
complicated CNDO method, is proven to have
embodied well-structured chemical bonds,
though it is not as successful for pi
bonds as it is for sigma bonds.
(for the whole file, see
log.staticCalculation.
The antibond occupancies are sometimes negative because of the
overlap transformation. Because the antibonds constructed by the
NBO method are not exactly the eigenstates of sub density matrices,
the numbers listed above are unimportant.)
Raising the antibond energies by 12 eV like shown in the above
lines (equivalent to lifting those in-gap electronic states
up to the band), we found that the energy gap becomes about
8 eV, compared with about 4 eV before changing. This is, apparently,
a significant improvement for the gap problem.
Major product of the EH
The EH/ET code more or less succeeds in producing a satisfactory
picture of bond structure for proteins, and, despite of the many
numerical approximations it has employed, succeeds as well in maintaining
the numerical stability and the conservation law of particle number.
The C-shell script that controls the job flow
This is the
script to create the input time series of structure
(in PDB) for EH/ET simulation. The script converts the DCD
segments and concatenates the segments into an integrated
time series.
#!/bin/csh -f
set name_of_protein = azur
set dir_of_strk = /opt/data/sdata/121-126/330ps
set dir_of_pdb = /home/xie/bet/strk
@ number_of_pdb_per_seg = 50
@ index_of_pdb = 1
@ iseg = 1
@ total = 1000
#check if the directory is correct
set a = `pwd`
if ( $a != $dir_of_pdb ) then
echo "We are not in the right starting directory!"
exit
endif
label_100:
echo $index_of_pdb, $iseg
cd $dir_of_strk
cp $name_of_protein.$iseg.tar.gz $dir_of_pdb
cd $dir_of_pdb
echo $name_of_protein.$iseg.tar.gz
gunzip $name_of_protein.$iseg.tar.gz
tar -xvf $name_of_protein.$iseg.tar
rm -f $name_of_protein.$iseg.tar
/bin/ls -1 *.pdb | sort -n -k 1.1 > tempo
@ j = ( $iseg - 1 ) * $number_of_pdb_per_seg
foreach ipdb ( `awk '{print $1}' tempo`)
@ j ++
if($j < 10) then
mv $ipdb $name_of_protein.0000$j
else if($j >= 10 && $j < 100) then
mv $ipdb $name_of_protein.000$j
else if($j >= 100 && $j < 1000) then
mv $ipdb $name_of_protein.00$j
else if($j >= 1000 && $j < 10000) then
mv $ipdb $name_of_protein.0$j
else if($j >= 10000 && $j < 100000) then
mv $ipdb $name_of_protein.$j
endif
echo $iseg, $ipdb, $j
end
@ iseg ++
if ( $iseg <= 60 ) goto label_100
exit
Because we have integrated the electronic structure calculation
part with the time-dependent electronic dynamics one, the script
is very simple:
#!/bin/csh -f set dir_of_main = /home/xie/bet #check if the directory is correct set a = `pwd` if ( $a != $dir_of_main ) then echo "We are not in the right starting directory!" exit endif @ irun = 0 label_100: bet_ns < restart_yes @ irun ++ if ( $irun < 150 ) goto label_100 exit
© 2000, Qian Xie