The following is the code snippet taken from the BONDO module.
This snippet does the job of CNDO's SCF interation. The lower case
and indented part of code is our modification. We have removed the
back transformation from BO to AO when z=it+1 from the
original code in
scfclo.f in order to reduce computational time. This back
transformation was for checking the NBO transformation. It is
not necessary for our calculation.
Important variables and arrays (Please refer to the CNDO formulism for details)
it: Number of iterations;
accel: To use Georgios's acceleration scheme or not;
energy: Total electronic energy;
q(nbszr): Coulomb integral;
g(nazr,nazr): Gamma matrix;
b(nbszr,nbszr): Density matrix;
a(nbszr,nbszr): Fock matrix;
cxie the density matrix input 'dmold' has its dimension written in
c the first place of the file, if this dimension identifier is
c different from the current dimension, the program will not take
c the input of 'dmold', instead it will use the default density
c matrix built using the Huckel model.
c a 'dmold' has to be present if the boolean flag 'accel' is
c set true, it doesn't matter if or not this 'dmold' is right.
c when 'accel' is set true, the program will ask for an input
c file 'dmold', if there is no such file, it will throw an error
c and stop.
c
if(accel) then
open(31,file='dmold')
read(31,*) n0
if(n0.ne.n) goto 1080
do idm = 1 , n
read(31,'(500f10.5)')(b(idm,jdm),jdm=1,n)
enddo
close(31)
endif
1080 continue
c GRAND SCF-ITERATIONS LOOP, Z.LT.it
c Z.LE.it SCF ITERATIONS BEFORE CONVERGENCE
c Z.EQ.it1 B.O. TRANSFORMATION, BASIS-SET TRUNCATION
c Z.EQ.it2 DEORTHOGONALIZATION AND EXIT
10 CONTINUE
IF(Z.NE.it) GO TO 20
WRITE(6,1600)
CALL EXIT
20 CONTINUE
Z = Z+1
if(z.eq.it2) goto 170
ENERGY = 0.D0
C
C CONSTRUCT FOCK MATRIX (IN A) FROM DENSITY MATRIX (IN B)0
C ...TRANSFER CORE HAMILTONIAN TO LOWER TRIANGLE OF A...
C
DO 30 I=1,N
A(I,I)=Q(I)
DO 30 J=I,N
30 A(J,I)=A(I,J)
DO 40 I=1,N
II=U(I)
A(I,I)=A(I,I)-B(I,I)*G(II,II)*0.5D0
DO 40 K=1,N
JJ=U(K)
40 A(I,I)=A(I,I)+B(K,K)*G(II,JJ)
NM=N-1
DO 50 I=1,NM
II=U(I)
LL=I+1
DO 50 J=LL,N
JJ=U(J)
50 A(J,I)=A(J,I)-B(J,I)*G(II,JJ)*0.5D0
C INDO MODIFICATION
IF (OPTION.EQ.CNDO) GO TO 100
write(6,*) ' INDO'
60 DO 90 II=1,NATOMS
K=AN(II)
I=LLIM(II)
IF (K.EQ.1) GO TO 90
70 PAA=B(I,I)+B(I+1,I+1)+B(I+2,I+2)+B(I+3,I+3)
A(I,I)=A(I,I)-(PAA-B(I,I))*G1(K)/6.D0
DO 80 J=1,3
A(I+J,I+J)=A(I+J,I+J)-B(I,I)*G1(K)/6.D0-(PAA-B(I,I))*7.D0*F2(K)/5
: 0.D0+B(I+J,I+J)*11.D0*F2(K)/50.D0
80 A(I+J,I)=A(I+J,I)+B(I,I+J)*G1(K)/2.D0
I1=I+1
I2=I+2
I3=I+3
A(I2,I1)=A(I2,I1)+B(I2,I1)*11.D0*F2(K)/50.D0
A(I3,I1)=A(I3,I1)+B(I3,I1)*11.D0*F2(K)/50.D0
A(I3,I2)=A(I3,I2)+B(I3,I2)*11.D0*F2(K)/50.D0
90 CONTINUE
C
C FOCK MATRIX COMPLETE; EVALUATE AND PRINT OUT ENERGY
C
100 CONTINUE
DO 110 I=1,N
110 ENERGY=ENERGY+0.5D0*B(I,I)*(A(I,I)+Q(I))
DO 120 I=1,NM
LL=I+1
DO 120 J=LL,N
120 ENERGY=ENERGY+B(I,J)*(A(I,J)+A(J,I))
IF((NLIST.EQ.0).AND.(Z.EQ.it2))GO TO 170
IF((Z.EQ.it2).AND.(IPR.GE.1)) WRITE(6,1200)
IF(IPR.GE.3) WRITE(6,1300) ENERGY
IF(DABS(ENERGY-OLDENG).GE.0.00005D0) GO TO 170
C
C ENERGY SATISFIED; SET Z=it1 FOR FINAL TWO PASSES
C
130 IF(Z.LE.IT)Z=it1
140 IF((Z.LT.it2).AND.(IPR.GE.1)) WRITE(6,1400)
IF(IPR.GE.1) WRITE(6,1300) ENERGY
IF(Z.EQ.it2)GO TO 170
C
C Z = it1 TRANSFORM, DELETE, DIAGONALIZE, PRINT B.O. EIGENVECTORS
C
IF((NLIST.EQ.0).AND.(IPR.LE.1)) GO TO 310
cxie save the convergent density matrix into a file 'dmold'
c for use in the next conformation. this was suggested by georgios,
c and it turns out to speed up BONDO calculation for a series of
c conformations taken from an MD trajectory
c
open(31,file='dmold')
rewind 31
write(31,*) n
do iof = 1 , n
write(31,'(500f10.5)')(b(iof,jof),jof=1,n)
enddo
close(31)
IF(IWLCAO.GT.0) GO TO 160
IF(IWNAT.EQ.0) GO TO 150
In our calculation the convergence criterion is set
as |E_new-E_old|<0.00005.
CNDO subroutines
We have written
makefile for HP-UX,
and makefile for Linux.
HP-UX has a native Fortran compiler. Linux has its peculiar
g77 complier. The makefiles were
written to serve compilation on different platforms. All the subroutines
for the CNDO/BONDO calculation can be found here:
HP-UX version and
Linux version. For CNDO calculation on the
atomic basis, just set the corresponding flag (see
input_lcao).
The Linux version fixes
a few bugs that only exist on the Linux platform (presumably owing
to the Linux Fortran compiler).
A list of some other BONDO subroutines and their functionalities
can be found here.
This script has been tested on HP-UX and Linux, however, there
is no warranty that it works as well on other UNIX.
In the case of BONDO calculation, the script extracts information
about the time series for the HOMO/LUMO( The script breaks the job
in the following three cases: (1) Directory error; (2)
BONDO module error (e.g. SCF energy not convergent within given
number of iterations); (3) NBO inappropriately constructed, namely,
the maximum antibond density exceeds 0.5.
© 2000, Qian Xie
Major CNDO product
Our code accepts only the PDB format for the structure input
file, and produces the following information about the electronic
structure
C-shell script that controls the job flow
Here we present a template script for obtaining the Hamiltonian
time series and extracting some other information
for a given segment of MD trajectory. The trajectory
segment is stored in UNIX tar and compression format. In order
to run this script, the excutable bondo should be
copied to the working directory (dir_of_et).
The user should also prepare two separate input files for the
BONDO module, input_lcao and input_lcbo,
which list the inputs for the
linear-combination-of-atomic-orbital and
the linear-combination-of-bond-orbital calculations, respectively.
#!/bin/csh -f
# compress_hamiltonian = 1 , compress ; = 0 , not compress (default)
# lcao = 1 , atomic orbital (default) ; = 0 , natural bond orbital
#
set initial_time = `date`
set name_of_protein = azur
set dir_of_et = /home/xie/ed8
set dir_of_bondo = /home/xie/ed8
set dir_of_strk = /home/xie/strk3
set dir_of_bondo_pdb = /home/xie/ed8/pdb
set dir_of_ham = /home/xie/ed8/hami
@ index_of_ladder = 1
@ lcao = 0
@ compress_hamiltonian = 0
@ number_of_pdb_per_seg = 50
@ index_of_seg = 0
/bin/date
/bin/pwd
/bin/rm -f temp1 temp2 temp3 temp4 temp5 tempo energy.tot
/bin/rm -f ./temp/temp*
# make the working directories if they do not exist.
# the directory to accomodate the hamiltonian matrices
if ( -e ./hami ) then
# do nothing at all.
else
mkdir hami
endif
# the directory which the pdb files are uncompressed to
if ( -e ./pdb ) then
# do nothing at all.
else
mkdir pdb
endif
# the directory where temporary files are stored
if ( -e ./temp ) then
# do nothing at all.
else
mkdir temp
endif
#check if the directory is correct
set a = `pwd`
if ( $a != $dir_of_et ) then
/bin/echo We are not in the right starting directory!
exit
endif
label_100:
@ iseg = 1
@ index_of_seg = $iseg
/bin/echo $iseg
# get pdb files
cd $dir_of_strk
cp $name_of_protein.$iseg.tar.gz $dir_of_bondo_pdb
cd $dir_of_bondo_pdb
rm -f *.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
# bondo/cndo calculation starts
cd $dir_of_bondo
if ( -e total ) /bin/rm total
if ( -e tempo ) /bin/rm tempo
cd $dir_of_bondo_pdb
/bin/ls -1 *.pdb | sort -n -k 1.1 > ../tempo
cd $dir_of_bondo
# check if the directory is correct
set a = `pwd`
if ( $a != $dir_of_bondo ) then
echo We are not in the BONDO directory!
exit
endif
if ( $lcao == 1 ) then
cp -f input_lcao input
else
cp -f input_lcbo input
endif
@ j = ( $iseg - 1 ) * $number_of_pdb_per_seg
foreach ipdb ( `awk '{print $1}' tempo`)
cp $dir_of_bondo_pdb/$ipdb proton.pdb
./bondo > log
# if BONDO calculation fails, stop anyway.
@ flag1 = `awk '{print $1}' flag_of_bondo`
if ( $flag1 == 0 ) then
echo BONDO error!
exit
endif
@ j ++
if($j < 10) then
cp fock fock.0000$j
if ( $compress_hamiltonian == 1 ) gzip fock.0000$j
else if($j >= 10 && $j < 100) then
cp fock fock.000$j
if ( $compress_hamiltonian == 1 ) gzip fock.000$j
else if($j >= 100 && $j < 1000) then
cp fock fock.00$j
if ( $compress_hamiltonian == 1 ) gzip fock.00$j
else if($j >= 1000 && $j < 10000) then
cp fock fock.0$j
if ( $compress_hamiltonian == 1 ) gzip fock.0$j
else if($j >= 10000 && $j < 100000) then
cp fock fock.$j
if ( $compress_hamiltonian == 1 ) gzip fock.$j
endif
mv fock.* $dir_of_ham
echo $index_of_seg, $ipdb, $j
grep -i 'total antibond density' log >> temp3
grep -i 'electron energy' log >> temp4
grep -i 'total energy' log >> energy.tot
grep -i 'homo' log >> temp1
grep -i 'lumo' log >> temp2
grep -i 'total antibond density' log > temp5
if ( $index_of_ladder == 1 ) then
grep -i '135. lp' log >> ./temp/temp_135
grep -i '110. bd' log >> ./temp/temp_110
endif
if ( $lcao == 1 ) goto label_200
set occmax_antibond = `awk '{print $11*10000}' temp5`
echo $occmax_antibond
if ( $occmax_antibond >= 5000 ) then
echo Maximum antibond occupancy too much!
exit
endif
label_200:
end
awk '{print NR,$4*2*13.6}' temp1 > ! homo.dat
awk '{print NR,$4*2*13.6}' temp2 > ! lumo.dat
awk '{print NR,$5,$11}' temp3 > ! antiden.tot
awk '{print NR,$3}' temp4 > ! elec.ene
cd $dir_of_bondo
label_010:
set final_time = `date`
echo "job started at $initial_time"
echo "job finished at $final_time"
exit
homo.dat, lumo.dat),
the total antibond density (antiden.tot), the total
electronic energy (elec.ene), and the detailed information
for given bonds (for example, in the above script, information about #135
lonepair and #110 bond were requested).