c$debug c-------------------------------------------------------------------------- c c EVAL - evaluation of bond lengths and bond angles from c list of cartesian coordinates with errors c c-------------------------------------------------------------------------- c Reedit any .EVA file: c c 1/ The first line with a numeric value in the first six columns c contains the number of Cartesians c 2/ Lines with Cartesians are fixed format, first six characters are c treated as a descriptor c 3/ The lines following Cartesians contain definitions of internals c in fixed format where: c - the first number defines 1=bond, 2=angle, 3=dihedral angle c - the following numbers refer to the sets of Cartesians c These lines are read up to the EOF c c C version 21.V.1999 ----- Zbigniew KISIEL ----- C C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- c parameter (maxat=100) implicit real*8 (a-h,o-z) character filnam*30,namat(maxat)*6 real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord c write(*,'(1x//)') WRITE(*,3344) 3344 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | E V A L - Internals and their errors from Cartesians' * ,T79,'|'/ * ' |',76(1H_),'|'/' version 21.V.1999',T64,'Zbigniew KISIEL'/) c 50 write(*,1)'File name containing data: ' 1 format(1x/1x,a,$) READ(*,'(A)',err=50)FILNAM c open(2,file=filnam,STATUS='OLD',err=50) c open(3,file='eval.out',status='unknown') write(3,3344) write(*,'(1x)') c c...read Cartesians c 15 read(2,'(i5)',err=15)natoms if(natoms.lt.1.or.natoms.gt.maxat)goto 15 do 16 j=1,natoms read(2,6,err=4)namat(j),(cord(j,k),ecord(j,k),k=1,3) 16 continue 6 format(a6,6f10.5) c write(3,'(1x/'' INPUT CARTESIANS:''/)') DO 11 J=1,NATOMS WRITE(*,10)NAMAT(J),(cord(j,k),ecord(j,k),k=1,3) WRITE(3,10)NAMAT(J),(cord(j,k),ecord(j,k),k=1,3) 11 continue 10 format(1x,a6,3(f10.5,f7.5)) c c...read definition of internal and do the appropriate calculation c write(*,'(1x)') write(3,'(1x/'' CALCULATED INTERNALS:''/)') c 7 read(2,8,err=4,end=4)ipar,ii,jj,kk,ll 8 format(5i3) if(ii.lt.1.or.ii.gt.natoms)goto 4 if(jj.lt.1.or.jj.gt.natoms)goto 4 if(ii.eq.jj)goto 4 c if(ipar.eq.1)then call bond(ii,jj,r,er) WRITE(*,9)' ',' ',namat(ii),namat(jj),R,ER WRITE(3,9)' ',' ',namat(ii),namat(jj),R,ER goto 7 endif 9 format(1x,4a6,' = ',F10.5,' +-',F8.5) c if(ipar.eq.2)then if(kk.lt.1.or.kk.gt.natoms)goto 4 if(ii.eq.kk.or.jj.eq.kk)goto 4 call angle(ii,jj,kk,an,ean) WRITE(*,9)' ',namat(ii),namat(jj),namat(kk),an,ean WRITE(3,9)' ',namat(ii),namat(jj),namat(kk),an,ean goto 7 endif c if(ipar.eq.3)then if(kk.lt.1.or.kk.gt.natoms)goto 4 if(ll.lt.1.or.ll.gt.natoms)goto 4 if(ii.eq.kk.or.ii.eq.ll)goto 4 if(jj.eq.kk.or.jj.eq.ll.or.kk.eq.ll)goto 4 call dangle(ii,jj,kk,ll,dan,edan) WRITE(*,9)namat(ii),namat(jj),namat(kk),namat(ll),dan WRITE(3,9)namat(ii),namat(jj),namat(kk),namat(ll),dan goto 7 endif C GOTO 7 C 4 close(3) c STOP END c c-------------------------------------------------------------------------- c subroutine bond(II,JJ,r,er) parameter (maxat=100) implicit real*8 (a-h,o-z) common /atoms/cord(maxat,3),ecord(maxat,3) a = cord(ii,1) da=ecord(ii,1) b = cord(ii,2) db=ecord(ii,2) c = cord(ii,3) dc=ecord(ii,3) a1 = cord(jj,1) da1=ecord(jj,1) b1 = cord(jj,2) db1=ecord(jj,2) c1 = cord(jj,3) dc1=ecord(jj,3) DELA=A1-A DELB=B1-B DELC=C1-C r=dsqrt( DELA**2 + DELB**2 + DELC**2 ) C EDELA=DSQRT(DA1**2+DA**2) EDELB=DSQRT(DB1**2+DB**2) EDELC=DSQRT(DC1**2+DC**2) C ER=DSQRT( (DELA**2*EDELA**2 + DELB**2*EDELB**2 * +DELC**2*EDELC**2)/R**2) c return end C c-------------------------------------------------------------------------- c subroutine ANGLE(ii,jj,kk,AN,EAN) c c The angle is that opposite r3 and is returned in degrees c parameter (maxat=100) implicit real*8 (a-h,o-z) common /atoms/cord(maxat,3),ecord(maxat,3) C call bond(ii,jj,r1,dr1) call bond(jj,kk,r2,dr2) call bond(ii,kk,r3,dr3) an=57.2958d0*dACOS((r1**2+r2**2-r3**2)/(2.d0*r1*r2)) c ean= dsqrt( $ -82070217441.d0*dr1**2*(r1**2-r2**2+r3**2)**2/(25000000.d0 $ *r1**2*(r1**4-2.d0*r1**2*(r2**2+r3**2)+r2**4-2.d0*r2**2*r3**2 $ +r3**4))-82070217441.d0*dr2**2*(r1**2-r2**2-r3**2)**2/ $ (25000000.d0*r2**2*(r1**4-2.d0*r1**2*(r2**2+r3**2)+r2**4- $ 2.d0*r2**2*r3**2+r3**4))-82070217441.d0*dr3**2*r3**2/ $ (6250000.d0*(r1**4-2.d0*r1**2*(r2**2+r3**2)+r2**4- $ 2.d0*r2**2*r3**2+r3**4)) ) c RETURN END C c-------------------------------------------------------------------------- c subroutine DANGLE(ii,jj,kk,ll,dan,edan) C C Dihedral angle (ii.jj.kk.ll), returned in degrees C parameter (maxat=100) implicit real*8 (a-h,o-z) real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord CONV=1.745329252D-2 C A1=cord(ii,1) A2=cord(ii,2) A3=cord(ii,3) B1=cord(jj,1) B2=cord(jj,2) B3=cord(jj,3) C1=cord(kk,1) C2=cord(kk,2) C3=cord(kk,3) D1=cord(ll,1) D2=cord(ll,2) D3=cord(ll,3) c c ANGL1=ANGLE(A1,A2,A3,B1,B2,B3,C1,C2,C3)-90.D0 c ANGL2=ANGLE(B1,B2,B3,C1,C2,C3,D1,D2,D3)-90.D0 c R12S=BONDL(A1,A2,A3,B1,B2,B3)*DSIN(ANGL1*CONV) c R23=BONDL(B1,B2,B3,C1,C2,C3) c R34S=BONDL(C1,C2,C3,D1,D2,D3)*DSIN(ANGL2*CONV) C call angle(ii,jj,kk,ANGL1,EANGL1) ANGL1=ANGL1-90.d0 call angle(jj,kk,ll,ANGL2,EANGL2) ANGL2=ANGL2-90.d0 call bond(ii,jj,R12,ER12) R12S=R12*DSIN(ANGL1*CONV) call bond(jj,kk,R23,er23) call bond(kk,ll,R34,er34) R34s=R34*DSIN(ANGL2*CONV) c P1=A1-B1+R12S*(C1-B1)/R23 P2=A2-B2+R12S*(C2-B2)/R23 P3=A3-B3+R12S*(C3-B3)/R23 Q1=D1-C1-R34S*(C1-B1)/R23 Q2=D2-C2-R34S*(C2-B2)/R23 Q3=D3-C3-R34S*(C3-B3)/R23 C BL1=DSQRT(P1**2+P2**2+P3**2) BL2=DSQRT(Q1**2+Q2**2+Q3**2) DOTP=P1*Q1+P2*Q2+P3*Q3 C DAN=DOTP/(BL1*BL2) IF(DABS(DAN).GT.1.D0)DAN=DSIGN(1.D0,DAN) DAN=DACOS(DAN)/CONV Edan=0.0d0 C RETURN END c c-------------------------------------------------------------------------- c--------------------------------------------------------------------------