Re: CCL:file format?
On Wed, 28 Feb 2001, Rajarshi Guha wrote:
> I was using the software package DTMM (for windows) to draw some
> molecules.
> Could anybody tell me what file format DTMM uses?
The answer appears to be 'a very fussy one'. My experience is that the
file format, while it is plain text, has to have the numbers in strictly
defined column positions.
For a modest sized molecule, it is easy to get from dtmm format by using
a text editor to hand edit the file to .xyz format, which can then be
translated using babel or molden. You just need a text editor which
preferably can do rectangular block copy, delete etc., i.e. is designed for
column-oriented text.
To go to dtmm format, I eventually became desperate enough to write a
program, todtmm.for, which I attach as plain text. This reads a coordinates
file which is in relatively free format, so can easily be generated from
an .xyz file by hand editing. This format is described in the program
comments.
My coordinates file format was originally designed for input to my
programs bangl and tangl, to tabulate bond lengths and angles, and
torsion angles, respectively. I also attach these programs, in case
they are of any use to you. They produce identical results to the
distance or angle tables in Gaussian output, or individual measurements
using molden (but different to Rasmol or the Chime plugin, which seem to
suffer from rounding errors).
No guarantees given! Bangl and tangl have been in use for many years,
but todtmm is much less tested, as I have now given up using dtmm in favour
of software which can read mdl .mol files etc.
Bruce Tattershall
Bruce.Tattershall "at@at" newcastle.ac.uk
Dr. B.W. Tattershall Department of Chemistry
University of Newcastle Newcastle upon Tyne
England
PROGRAM TODTMM
C=== READS .CRD FILE, SUITABLE FOR BANGL;
C=== WRITES .MOL FILE, SUITABLE FOR DTMM
C===
C=== PROGRAM CREATED BY MODIFYING BANGL
C===
C=== INPUT FILE (STANDARD INPUT) CONSISTS OF:
C=== TITLE (A79)
C=== NO. OF ATOMS
C=== MAX. LENGTH OF BONDS
C=== THEN NATOMS PAIRS OF LINES:
C=== ATOM LABEL (A8)
C=== X Y Z (FREE FORMAT)
C===
C=== HANDLES UP TO 100 ATOMS, 400 BONDS, COORD. NO. 10
C===
REAL X(100),Y(100),Z(100)
DOUBLE PRECISION BLEN(400),BLEN2(400),DIST2,
& XI,YI,ZI,BLJ,BL2J,BONDMX,BONDM2
INTEGER IBNDED(400),JBNDED(400),LIGND(10)
CHARACTER SYMBOL(100)*8,SYMBI*8
C===
DIST2(JATOM)=(XI-X(JATOM))**2+(YI-Y(JATOM))**2+(ZI-Z(JATOM))**2
C===
MATOMS=400
MBONDS=100
MLIGS=10
C===
ONE=1.0
RNINTY=90.0
C===
C=== READ IN DATA
READ(*,10)
10 FORMAT(A79)
READ(*,30)NATOMS
30 FORMAT(I3)
IF(NATOMS.LE.MATOMS)GOTO 50
WRITE(*,40)MATOMS
40 FORMAT(/' TOO MANY ATOMS, LIMIT IS ',I3)
STOP
50 READ(*,60)BONDMX
60 FORMAT(F10.0)
DO 80 IATOM=1,NATOMS
READ(*,70,END=240)SYMBOL(IATOM)
70 FORMAT(A8)
80 READ(*,*,END=240,ERR=240)X(IATOM),Y(IATOM),Z(IATOM)
C===
NBONDS=0
BONDM2=BONDMX**2
C=== LOOP TO STORE BOND LENGTHS BLEN(), THEIR SQUARES BLEN2(),
C=== AND THE INDICES OF ATOMS THEY ARE BETWEEN, IBNDED() AND JBNDED()
DO 110 IATOM=1,NATOMS-1
XI=X(IATOM)
YI=Y(IATOM)
ZI=Z(IATOM)
DO 110 JATOM=IATOM+1,NATOMS
BL2J=DIST2(JATOM)
IF(BL2J.GT.BONDM2)GOTO 110
IF(NBONDS.LT.MBONDS)GOTO 100
WRITE(*,90)MBONDS,BONDMX
90 FORMAT(/'TOO MANY BONDS; LIMIT IS ',I4/
& ' MAX. BOND LENGTH INPUT WAS',F7.3)
STOP
100 NBONDS=NBONDS+1
IBNDED(NBONDS)=IATOM
JBNDED(NBONDS)=JATOM
110 CONTINUE
C=== WRITE HEADER LINES
WRITE(*,112)ONE,ONE,ONE,RNINTY,RNINTY,RNINTY,NATOMS
112 FORMAT(T42,F5.3,T50,F5.3,T58,F5.3/T24,F6.3,T32,F6.3,T40,F6.3/
& I4/' ')
C===
C=== FOR EACH ATOM, MAKE TEMPORARY ARRAYS
C=== INBNDS() OF INDEXES TO BONDS ARRAYS
C=== LIGND() OF INDEXES TO ATOMS ARRAYS
DO 230 IATOM=1,NATOMS
SYMBI=SYMBOL(IATOM)
DO 120 KBOND=1,NBONDS
IF(IBNDED(KBOND).EQ.IATOM)GOTO 130
120 CONTINUE
130 NLIGS=0
IF(KBOND.EQ.1)GOTO 145
C=== BONDS FROM EARLIER ATOMS WILL BE SPREAD OUT UP TO HERE
DO 140 IBOND=1,KBOND-1
IF(JBNDED(IBOND).NE.IATOM)GOTO 140
NLIGS=NLIGS+1
IF(NLIGS.GT.MLIGS)GOTO 160
LIGND(NLIGS)=IBNDED(IBOND)
140 CONTINUE
C===
C=== BONDS TO LATER ATOMS WILL BE BUNCHED HERE
145 DO 150 IBOND=KBOND,NBONDS
IF(IBNDED(IBOND).NE.IATOM)GOTO 180
NLIGS=NLIGS+1
IF(NLIGS.GT.MLIGS)GOTO 160
LIGND(NLIGS)=JBNDED(IBOND)
150 CONTINUE
GOTO 180
C===
160 WRITE(*,170)SYMBI,MLIGS,BONDMX
170 FORMAT(/' ATOM ',A8,': TOO MANY LIGANDS; LIMIT IS ',I3/
& ' MAX. BOND LENGTH INPUT WAS',F7.3)
STOP
C===
C=== WRITE ATOM RECORD
180 WRITE(*,190)IATOM,SYMBI,X(IATOM),Y(IATOM),Z(IATOM),
& (LIGND(ILIGND),ILIGND=1,NLIGS)
190 FORMAT(I4,1X,A2,3X,3F10.5,1X,4I4)
C===
230 CONTINUE
STOP
C===
240 WRITE(*,250)
250 FORMAT(/' BAD DATA')
STOP
C===
END
L PROGRAM BANGL
C=== CALCULATES BOND LENGTHS AND ANGLES FROM CARTESIAN COORDS.
C===
C=== INPUT FILE (STANDARD INPUT) CONSISTS OF:
C=== TITLE (A79)
C=== NO. OF ATOMS
C=== MAX. LENGTH OF BONDS
C=== THEN NATOMS PAIRS OF LINES:
C=== ATOM LABEL (A8)
C=== X Y Z (FREE FORMAT)
C===
C=== HANDLES UP TO 100 ATOMS, 400 BONDS, COORD. NO. 10
C===
REAL X(100),Y(100),Z(100)
DOUBLE PRECISION BLEN(400),BLEN2(400),DIST2,
& XI,YI,ZI,BLJ,BL2J,BONDMX,BONDM2
INTEGER IBNDED(400),JBNDED(400),LIGND(10),INBNDS(10)
CHARACTER SYMBOL(100)*8,SYMBI*8,SYMBJ*8,TITLE*79
C===
DIST2(JATOM)=(XI-X(JATOM))**2+(YI-Y(JATOM))**2+(ZI-Z(JATOM))**2
C===
MATOMS=400
MBONDS=100
MLIGS=10
RTOD=180.0/ACOS(-1.0)
C===
C=== READ IN DATA
READ(*,10)TITLE
10 FORMAT(A79)
WRITE(*,20)TITLE
20 FORMAT(1X,A79)
READ(*,30)NATOMS
30 FORMAT(I3)
IF(NATOMS.LE.MATOMS)GOTO 50
WRITE(*,40)MATOMS
40 FORMAT(/' TOO MANY ATOMS, LIMIT IS ',I3)
STOP
50 READ(*,60)BONDMX
60 FORMAT(F10.0)
DO 80 IATOM=1,NATOMS
READ(*,70,END=240)SYMBOL(IATOM)
70 FORMAT(A8)
80 READ(*,*,END=240,ERR=240)X(IATOM),Y(IATOM),Z(IATOM)
C===
NBONDS=0
BONDM2=BONDMX**2
C=== LOOP TO STORE BOND LENGTHS BLEN(), THEIR SQUARES BLEN2(),
C=== AND THE INDICES OF ATOMS THEY ARE BETWEEN, IBNDED() AND JBNDED()
DO 110 IATOM=1,NATOMS-1
XI=X(IATOM)
YI=Y(IATOM)
ZI=Z(IATOM)
DO 110 JATOM=IATOM+1,NATOMS
BL2J=DIST2(JATOM)
IF(BL2J.GT.BONDM2)GOTO 110
IF(NBONDS.LT.MBONDS)GOTO 100
WRITE(*,90)MBONDS,BONDMX
90 FORMAT(/'TOO MANY BONDS; LIMIT IS ',I4/
& ' MAX. BOND LENGTH INPUT WAS',F7.3)
STOP
100 NBONDS=NBONDS+1
IBNDED(NBONDS)=IATOM
JBNDED(NBONDS)=JATOM
BLEN2(NBONDS)=BL2J
BLEN(NBONDS)=SQRT(BL2J)
110 CONTINUE
C===
C=== FOR EACH ATOM, MAKE TEMPORARY ARRAYS
C=== INBNDS() OF INDEXES TO BONDS ARRAYS
C=== LIGND() OF INDEXES TO ATOMS ARRAYS
DO 230 IATOM=1,NATOMS
SYMBI=SYMBOL(IATOM)
DO 120 KBOND=1,NBONDS
IF(IBNDED(KBOND).EQ.IATOM)GOTO 130
120 CONTINUE
130 NLIGS=0
IF(KBOND.EQ.1)GOTO 145
C=== BONDS FROM EARLIER ATOMS WILL BE SPREAD OUT UP TO HERE
DO 140 IBOND=1,KBOND-1
IF(JBNDED(IBOND).NE.IATOM)GOTO 140
NLIGS=NLIGS+1
IF(NLIGS.GT.MLIGS)GOTO 160
INBNDS(NLIGS)=IBOND
LIGND(NLIGS)=IBNDED(IBOND)
140 CONTINUE
C===
C=== BONDS TO LATER ATOMS WILL BE BUNCHED HERE
145 DO 150 IBOND=KBOND,NBONDS
IF(IBNDED(IBOND).NE.IATOM)GOTO 180
NLIGS=NLIGS+1
IF(NLIGS.GT.MLIGS)GOTO 160
INBNDS(NLIGS)=IBOND
LIGND(NLIGS)=JBNDED(IBOND)
150 CONTINUE
GOTO 180
C===
160 WRITE(*,170)SYMBI,MLIGS,BONDMX
170 FORMAT(/' ATOM ',A8,': TOO MANY LIGANDS; LIMIT IS ',I3/
& ' MAX. BOND LENGTH INPUT WAS',F7.3)
STOP
C===
180 IF(NLIGS.EQ.0)GOTO 230
C=== OUTPUT BOND LENGTHS
WRITE(*,190)SYMBI,NLIGS
190 FORMAT(/' ATOM ',A8,': COORDINATION NUMBER ',I2/
& ' BOND LENGTHS:')
DO 200 ILIGND=1,NLIGS
JATOM=LIGND(ILIGND)
IBOND=INBNDS(ILIGND)
200 WRITE(*,210)SYMBI,SYMBOL(JATOM),BLEN(IBOND)
210 FORMAT(1X,A8,' -- ',A8,2X,F6.4)
C===
IF(NLIGS.EQ.1)GOTO 230
C=== OUTPUT BOND ANGLES
WRITE(*,220)SYMBI
220 FORMAT(' ANGLES ABOUT ATOM ',A8,':')
DO 230 JLIGND=1,NLIGS-1
JATOM=LIGND(JLIGND)
JBOND=INBNDS(JLIGND)
XI=X(JATOM)
YI=Y(JATOM)
ZI=Z(JATOM)
SYMBJ=SYMBOL(JATOM)
BL2J=BLEN2(JBOND)
BLJ=BLEN(JBOND)
DO 230 KLIGND=JLIGND+1,NLIGS
KATOM=LIGND(KLIGND)
KBOND=INBNDS(KLIGND)
ANGLE=RTOD*ABS(ACOS((BL2J+BLEN2(KBOND)-DIST2(KATOM))/
& (2*BLJ*BLEN(KBOND))))
WRITE(*,229)SYMBJ,SYMBI,SYMBOL(KATOM),ANGLE
229 FORMAT(1X,A8,' -- ',A8,' -- ',A8,2X,F7.3)
C===
230 CONTINUE
STOP
C===
240 WRITE(*,250)
250 FORMAT(/' BAD DATA')
STOP
C===
END
PROGRAM TANGL
C=== CALCULATES TORSION ANGLES A-B-C-D FROM CARTESIAN COORDS.
C===
C=== INPUT FILE (STANDARD INPUT) IS THE SAME AS FOR BANGL, AND
C=== CONSISTS OF:
C=== TITLE (A79)
C=== NO. OF ATOMS
C=== MAX. LENGTH OF BONDS
C=== THEN NATOMS PAIRS OF LINES:
C=== ATOM LABEL (A8)
C=== X Y Z (FREE FORMAT)
C===
C=== HANDLES UP TO 100 ATOMS, 400 BONDS, COORD. NO. 10
C===
C***
COMMON X,Y,Z,NBONDS,IBNDED,JBNDED,MLIGS,BONDMX
DOUBLE PRECISION BONDMX
REAL X(100),Y(100),Z(100)
INTEGER IBNDED(400),JBNDED(400)
C***
DOUBLE PRECISION RTOD,DIST2,XI,YI,ZI,BL2J,BONDM2,ANGLA,ANGLT,
& ANGLD,VECBA(3),VECBC(3),VECCB(3),VECCD(3),VECCR1(3),VECCR2(3),
& VECCR3(3),VDOT,VANGL
INTEGER LIGNDB(10),LIGNDC(10)
CHARACTER SYMBOL(100)*8,SYMBA*8,SYMBB*8,SYMBC*8,SYMBD*8,TITLE*79
C===
DIST2(JATOM)=(XI-X(JATOM))**2+(YI-Y(JATOM))**2+(ZI-Z(JATOM))**2
C===
MATOMS=400
MBONDS=100
MLIGS=10
RTOD=DBLE(180.0)/DACOS(DBLE(-1.0))
C===
C=== READ IN DATA
READ(*,10)TITLE
10 FORMAT(A79)
WRITE(*,20)TITLE
20 FORMAT(1X,A79)
READ(*,30)NATOMS
30 FORMAT(I3)
IF(NATOMS.LE.MATOMS)GOTO 50
WRITE(*,40)MATOMS
40 FORMAT(/' TOO MANY ATOMS, LIMIT IS ',I3)
STOP
50 READ(*,60)BONDMX
60 FORMAT(F10.0)
DO 80 IATOM=1,NATOMS
READ(*,70,END=240)SYMBOL(IATOM)
70 FORMAT(A8)
80 READ(*,*,END=240,ERR=240)X(IATOM),Y(IATOM),Z(IATOM)
C===
NBONDS=0
BONDM2=BONDMX**2
C=== LOOP TO STORE THE INDICES OF ATOMS BONDS ARE BETWEEN,
C=== IBNDED() AND JBNDED()
DO 110 IATOM=1,NATOMS-1
XI=X(IATOM)
YI=Y(IATOM)
ZI=Z(IATOM)
DO 110 JATOM=IATOM+1,NATOMS
BL2J=DIST2(JATOM)
IF(BL2J.GT.BONDM2)GOTO 110
IF(NBONDS.LT.MBONDS)GOTO 100
WRITE(*,90)MBONDS,BONDMX
90 FORMAT(/'TOO MANY BONDS; LIMIT IS ',I4/
& ' MAX. BOND LENGTH INPUT WAS',F7.3)
STOP
100 NBONDS=NBONDS+1
IBNDED(NBONDS)=IATOM
JBNDED(NBONDS)=JATOM
110 CONTINUE
C===
C=== HEADING FOR TABLE
WRITE(*,115)
115 FORMAT(/,3X,'A',11X,'B',11X,'C',11X,'D',6X,'Layback A',2X,
& 'Torsion',2X,'Layback D')
C=== FOR EACH ATOM B
DO 230 IATOMB=1,NATOMS-1
SYMBB=SYMBOL(IATOMB)
C=== GET THE LIGAND SET FOR B
CALL LIGNDS(IATOMB,NLIGSB,LIGNDB)
IF(NLIGSB.LT.2)GOTO 230
C===
C=== FOR EACH ATOM C
DO 230 ILC=1,NLIGSB
IATOMC=LIGNDB(ILC)
C=== EXCLUDE BC BONDS DONE ALREADY
IF(IATOMC.LT.IATOMB)GOTO 230
SYMBC=SYMBOL(IATOMC)
C=== SPACE THE OUTPUT
WRITE(*,120)
120 FORMAT(' ')
C=== GET LIGAND SET FOR C
CALL LIGNDS(IATOMC,NLIGSC,LIGNDC)
C=== VECTOR B->C
CALL VBTWEN(IATOMB,IATOMC,VECBC)
C===
C=== FOR EACH ATOM A
DO 230 ILA=1,NLIGSB
C=== EXCLUDE ATOM C
IF(ILA.EQ.ILC)GOTO 230
IATOMA=LIGNDB(ILA)
SYMBA=SYMBOL(IATOMA)
C=== VECTOR B->A
CALL VBTWEN(IATOMB,IATOMA,VECBA)
C=== CALC CROSS PRODUCT BC x BA
CALL VCROSS(VECBC,VECBA,VECCR1)
C=== CALC LAYBACK ANGLE ABC-90
ANGLA=VANGL(VECBC,VECBA,RTOD)-9.D1
C===
C=== FOR EACH ATOM D
DO 230 ILD=1,NLIGSC
IATOMD=LIGNDC(ILD)
C=== EXCLUDE ATOM B
IF(IATOMD.EQ.IATOMB)GOTO 230
SYMBD=SYMBOL(IATOMD)
C=== VECTOR C->D
CALL VBTWEN(IATOMC,IATOMD,VECCD)
C=== CALC CROSS PRODUCT CD x CB
C=== REVERSE VECTOR BC
DO 150 I=1,3
150 VECCB(I)=-VECBC(I)
CALL VCROSS(VECCD,VECCB,VECCR2)
C=== CALC LAYBACK ANGLE BCD-90
ANGLD=VANGL(VECCD,VECCB,RTOD)-9.D1
C=== THE CROSS PRODUCTS NOW IN VECCR1 AND VECCR2 ARE IN THE PLANE
C=== PERPENDICULAR TO B->C AND ARE AT 90 DEG TO THE PROJECTIONS
C=== OF B->A AND C->D RESPECTIVELY ONTO THIS PLANE
C===
C=== NOW FIND DOT PRODUCT OF THESE AS A MEANS OF FINDING THE
C=== ANGLE BETWEEN THEM AND
C=== COMPARE DIRECTION OF CROSS PRODUCT OF THEM, WHICH WILL BE
C=== PARALLEL TO B->C, WITH B->C, TO GET SIGN OF ANGLE
C===
ANGLT=VANGL(VECCR1,VECCR2,RTOD)
IF(ANGLT.EQ.3.6D2)GOTO 210
CALL VCROSS(VECCR1,VECCR2,VECCR3)
ANGLT=DSIGN(ANGLT,VDOT(VECCR3,VECBC))
C===
WRITE(*,200)SYMBA,SYMBB,SYMBC,SYMBD,ANGLA,ANGLT,ANGLD
200 FORMAT(1X,A8,3(' -- ',A8),3(2X,F8.3))
GOTO 230
210 WRITE(*,220)SYMBA,SYMBB,SYMBC,SYMBD,ANGLA,ANGLD
220 FORMAT(1X,A8,3(' -- ',A8),2X,F8.3,10X,F8.3)
C===
230 CONTINUE
C===
STOP
C===
240 WRITE(*,250)
250 FORMAT(/' BAD DATA')
STOP
C===
END
C*****
SUBROUTINE LIGNDS(IATOM,NLIGS,LIGND)
C=== COUNT LIGANDS AND MAKE TEMPORARY ARRAY
C=== LIGND() OF INDEXES TO ATOMS ARRAYS
C***
COMMON X,Y,Z,NBONDS,IBNDED,JBNDED,MLIGS,BONDMX
DOUBLE PRECISION BONDMX
REAL X(100),Y(100),Z(100)
INTEGER IBNDED(400),JBNDED(400)
C***
INTEGER LIGND(10)
C===
DO 120 KBOND=1,NBONDS
IF(IBNDED(KBOND).EQ.IATOM)GOTO 130
120 CONTINUE
130 NLIGS=0
IF(KBOND.EQ.1)GOTO 145
C=== BONDS FROM EARLIER ATOMS WILL BE SPREAD OUT UP TO HERE
DO 140 IBOND=1,KBOND-1
IF(JBNDED(IBOND).NE.IATOM)GOTO 140
NLIGS=NLIGS+1
IF(NLIGS.GT.MLIGS)GOTO 160
LIGND(NLIGS)=IBNDED(IBOND)
140 CONTINUE
C===
C=== BONDS TO LATER ATOMS WILL BE BUNCHED HERE
145 DO 150 IBOND=KBOND,NBONDS
IF(IBNDED(IBOND).NE.IATOM)GOTO 155
NLIGS=NLIGS+1
IF(NLIGS.GT.MLIGS)GOTO 160
LIGND(NLIGS)=JBNDED(IBOND)
150 CONTINUE
155 RETURN
C===
160 WRITE(*,170)SYMBI,MLIGS,BONDMX
170 FORMAT(/' ATOM ',A8,': TOO MANY LIGANDS; LIMIT IS ',I3/
& ' MAX. BOND LENGTH INPUT WAS',F7.3)
STOP
END
C*****
SUBROUTINE VBTWEN(IATOM,JATOM,VEC)
C=== LOOKS UP ATOM COORDS AND CALCULATES VECTOR BETWEEN THEM
C***
COMMON X,Y,Z,NBONDS,IBNDED,JBNDED,MLIGS,BONDMX
DOUBLE PRECISION BONDMX
REAL X(100),Y(100),Z(100)
INTEGER IBNDED(400),JBNDED(400)
C***
DOUBLE PRECISION VEC(3)
C===
VEC(1)=X(JATOM)-X(IATOM)
VEC(2)=Y(JATOM)-Y(IATOM)
VEC(3)=Z(JATOM)-Z(IATOM)
RETURN
END
C*****
SUBROUTINE VCROSS(VEC1,VEC2,VEC3)
C=== CALCULATES THE CROSS PRODUCT
C=== VEC3 = VEC1 x VEC2
C===
DOUBLE PRECISION VEC1(3),VEC2(3),VEC3(3)
C===
VEC3(1)=VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
VEC3(2)=VEC1(3)*VEC2(1)-VEC1(1)*VEC2(3)
VEC3(3)=VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
RETURN
END
C*****
DOUBLE PRECISION FUNCTION VANGL(VEC1,VEC2,RTOD)
C=== CALCULATES ANGLE IN DEGREES BETWEEN TWO VECTORS
C=== RTOD IS THE RADIANS TO DEGREES CONVERSION FACTOR
DOUBLE PRECISION VEC1(3),VEC2(3),VDOT,VMAG,RTOD,COSANG,VLEN1,VLEN2
VLEN1=VMAG(VEC1)
VLEN2=VMAG(VEC2)
VANGL=3.6D2
IF(VLEN1.LT.1.D-12.OR.VLEN2.LT.1.D-12)GOTO 10
COSANG=VDOT(VEC1,VEC2)/(VLEN1*VLEN2)
IF(COSANG.GT.1.D0)COSANG=1.D0
IF(COSANG.LT.-1.D0)COSANG=-1.D0
VANGL=RTOD*DACOS(COSANG)
10 RETURN
END
C*****
DOUBLE PRECISION FUNCTION VDOT(VEC1,VEC2)
C=== CALCULATES THE DOT PRODUCT VEC1.VEC2
DOUBLE PRECISION VEC1(3),VEC2(3)
VDOT=0.0
DO 10 I=1,3
10 VDOT=VDOT+VEC1(I)*VEC2(I)
RETURN
END
C*****
DOUBLE PRECISION FUNCTION VMAG(VEC)
C=== CALCULATES MAGNITUDE OF VECTOR VMAG
DOUBLE PRECISION VEC(3),DZERO,COMP,SUM
DZERO=DBLE(0.0)
SUM=DZERO
DO 10 I=1,3
COMP=VEC(I)
IF(DABS(COMP).LT.1.D-12)GOTO 10
SUM=SUM+COMP*COMP
10 CONTINUE
VMAG=DZERO
IF(SUM.LE.DZERO)GOTO 20
VMAG=DSQRT(SUM)
20 RETURN
END
C*****
E