c-------------------------------------------------------------------------- c c EVAL - evaluation of distances and angles from c a list of Cartesian coordinates with errors c c-------------------------------------------------------------------------- C c Usage: C c Reedit any previous .EVA file. Note that lines that begin with a c recognised comment flag (#,!,%) are ignored so that c comment lines can be inserted anywhere in the data file C C C The significant lines in the data are to be in the following C sequence and should contain: c c 1/ The number of lines that contain the (atomic) Cartesians C C 2/ The lines containing the Cartesians in the form: C DESCR,X,dx,Y,dY,Z,dZ C where: DESCR = six character alphanumeric descriptor C X,dX = the Cartesian and its error etc. C C 3/ The lines with definitions of coordinates to be evaluated C from the (atomic) Cartesians. These can be the standard structural C internal coordinates or some other geometrical parameters. C The possibilities are: C C 1,i,j C 2,i,j,k C 3,i,j,k,l C C -1,i,j,k,l C -2,i,j,k,l,m,n C -3,i,j,k,l,m C -4,i,j,k,l C -5,i,j,k,l C -6,i,j,k C C where: C C The first number is a parameter identifier specifying: C 1 = bond, 2 = angle, 3 = dihedral angle C -1 = plane containing centres i,j,k C -2 = angle between planes i,j,k and l,m,n C -3 = angle between plane i,j,k and line l,m C -4 = angle between line i,j and line k,l C -5 = the (perpendicular) distance between plane i,j,k and point l C -6 = the (perpendicular) distance between line j,j and point l. C c The numbers following the parameter identifier refer to the declared c Cartesians which are numbered in the order of declaration. c c The line defining a parameter to de evaluated can also conclude with c a comment, which has to begin with the '!' character. This comment will c be echoed to output. c c The parameter definition lines are read up to the end of file or to c a line beginning with the string: end c c c WARNING: c c The EVAL uncertainties in internals are evaluated without access to c the full variance-covariance matrix from the fit used to evaluate c the cited uncertainties in the input Cartesians. Unit diagonal c correlation matrix is assumed. c c Thus the uncertainties from EVAL may differ significantly (but typically c by not more than by 30% either way) from uncertainties in explicitly fitted c internals corresponding to the input Cartesians. c c Uncertainties in evaluation of the equationof a plane may differ depending c on the order of declared points. The first point is used twice so it is c recommended to declare the most precise point as point i. c c-------------------------------------------------------------------------- c Sample input file (first three characters in each line added for clarity) c c ! c ! pyrimidine, r_s coordinates c ! c 6 c N' -1.1982 0.0013 0.6848 0.0022 0.0 0.0 c C(2) 0.0 0.0 1.2791 0.0012 0.0 0.0 c N 1.1982 0.0013 0.6848 0.0022 0.0 0.0 c C(4) 1.1843 0.0013 -0.6468 0.0023 0.0 0.0 c C(5) 0.0 0.0 -1.3804 0.0011 0.0 0.0 c C(6) -1.1843 0.0013 -0.6468 0.0023 0.0 0.0 c ! c ! Bond lengths c ! c 1 2 3 c 1 3 4 c 1 4 5 c ! c ! Bond angles c ! c 2 1 2 3 c 2 2 3 4 ! comments can also be placed here c 2 3 4 5 c 2 4 5 6 c end c c C version 16.V.2020 ----- 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 c Modification history: c c ca 95: Created ad hoc in order to process results from KRA, c several minor mods made since then C 22.10.12: More robust input, enhanced commenting, unified error derivation C including the much demanded error in an evaluated dihedral c 19.06.17: Evaluation of additional parameters specified by negative indices c 16.05.20: Removing code inconsistencies revealed by Fortran2018 warnings C_____________________________________________________________________________ c c parameter (maxat=500) implicit real*8 (a-h,o-z) character filnam*100,namat(maxat)*6,errmes*50,line*100,strbuf*100, * filout*100 real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord CONV=1.7453292519943D-2 2 pi/360 c nlinp=0 write(*,'(1x//)') WRITE(*,3344) 3344 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | E V A L - Internals and their errors from Cartesians' * ,T79,'|'/ * ' |',76(1H_),'|'/' version 20.V.2020',T64,'Zbigniew KISIEL'/) c write(*,3346) 3346 format(/' Available evaluations:'//5x, * 'D= 1,i,j P = -1,i,j,k,l A(L,L)= -4,i,j,k,l'/ * 5x, * 'A= 2,i,j,k A(P,P)= -2,i,j,k,l,m,n D(P,p)= -5,i,j,k,l'/ * 5x, * 'T= 3,i,j,k,l A(P,L)= -3,i,j,k,l,m D(L,p)= -6,i,j,k'/) 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',form='FORMATTED',err=50) filout=filnam(1:len_trim(filnam))//'.out' c open(3,file=filout(1:len_trim(filout)),status='unknown') write(3,3344) write(3,3345) 3345 format(1x/' WARNING:'// * 10x,'The EVAL uncertainties are evaluated by assuming that'/ * 10x,'the correlation matrix is a unit matrix. '// * 10x,'The EVAL uncertainties may thus differ significantly '/ * 10x,'(but typically by not more than 30% either way) '/ * 10x,'from uncertainties in explicitly fitted internals '/ * 10x,'corresponding to the input Cartesians. '//) write(*,'(1x)') c c...read the Cartesians c errmes='Problem with the number of Cartesian lines' 15 nlinp=nlinp+1 call test(2,inperr) if(inperr.eq.1)goto 15 read(2,*,err=2000)natoms if(natoms.lt.1.or.natoms.gt.maxat)goto 2000 c errmes='Problem reading a line of Cartesians' do 16 j=1,natoms 20 nlinp=nlinp+1 call test(2,inperr) if(inperr.eq.1)goto 20 read(2,'(a)')line read(line,'(a)',err=2000,end=2000)namat(j) line=line(7:) read(line,*,err=2000,end=2000)(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(f12.5,f9.5)) C c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...read a definition of an internal and perform the appropriate calculation c write(*,'(1x)') write(3,'(1x/'' CALCULATED INTERNALS:''/)') c errmes='Problem reading the declaration of an internal' C 7 nlinp=nlinp+1 call test(2,inperr) if(inperr.eq.1)goto 7 c read(2,'(a)',err=2000,end=4)line if(line(1:3).eq.'end')goto 4 c if(len_trim(line).le.3)goto 7 read(line,*,err=2000,end=2000)ipar if(ipar.ne.-1.and.ipar.ne.-2.and.ipar.ne.-3.and.ipar.ne.-4. * .and.ipar.ne.-5.and.ipar.ne.-6. * .and.ipar.ne.1.and.ipar.ne.2.and.ipar.ne.3)then write(*,'(1x/'' ERROR: illegal internal index = '',i3//)')ipar stop endif c if(ipar.eq.1)read(line,*,err=2000,end=2000)ipar,ii,jj distance if(ipar.eq.2)read(line,*,err=2000,end=2000)ipar,ii,jj,kk angle if(ipar.eq.3)read(line,*,err=2000,end=2000)ipar,ii,jj,kk,ll dihedral if(ipar.eq.-1)read(line,*,err=2000,end=2000)ipar,ii,jj,kk plane if(ipar.eq.-2)read(line,*,err=2000,end=2000) * ipar,ii,jj,kk,ll,mm,nn plane.plane angle if(ipar.eq.-3)read(line,*,err=2000,end=2000)ipar,ii,jj,kk,ll,mm plane.line angle if(ipar.eq.-4)read(line,*,err=2000,end=2000)ipar,ii,jj,kk,ll line.line angle if(ipar.eq.-5)read(line,*,err=2000,end=2000)ipar,ii,jj,kk,ll plane.point distance if(ipar.eq.-6)read(line,*,err=2000,end=2000)ipar,ii,jj,kk line.point distance 8 format(5i3) c errmes='Atom number out of range' if(ii.lt.1.or.ii.gt.natoms)goto 2000 if(jj.lt.1.or.jj.gt.natoms)goto 2000 if(ii.eq.jj)goto 2000 c iend=0 do 136 n=len_trim(line),5,-1 if(line(n:n).eq.'!')then iend=n goto 137 endif 136 continue 137 if(iend.eq.0)then strbuf='' else strbuf=line(iend:len_trim(line)) endif c c c------------------------------------------------------------------------------- c Distance c------------------------------------------------------------------------------- c if(ipar.eq.1)then c call bond(ii,jj,r) <---- call ebond(ii,jj,r,er) <---- c WRITE(*,9)' ',' ',namat(ii),namat(jj),R,ER, * strbuf(1:len_trim(strbuf)) WRITE(3,9)' ',' ',namat(ii),namat(jj),R,ER, * strbuf(1:len_trim(strbuf)) goto 7 endif 9 format(1x,4a6,' = ',F10.5,' +-',F8.5,5x,a) c c c------------------------------------------------------------------------------- c Angle c------------------------------------------------------------------------------- c if(ipar.eq.2)then if(kk.lt.1.or.kk.gt.natoms)goto 2000 if(ii.eq.kk.or.jj.eq.kk)goto 2000 c call angle(ii,jj,kk,an) <---- call eangle(ii,jj,kk,an,ean) <---- c WRITE(*,9)' ',namat(ii),namat(jj),namat(kk),an,ean, * strbuf(1:len_trim(strbuf)) WRITE(3,9)' ',namat(ii),namat(jj),namat(kk),an,ean, * strbuf(1:len_trim(strbuf)) goto 7 endif c c c------------------------------------------------------------------------------- c Dihedral angle c------------------------------------------------------------------------------- c if(ipar.eq.3)then if(kk.lt.1.or.kk.gt.natoms)goto 2000 if(ll.lt.1.or.ll.gt.natoms)goto 2000 if(ii.eq.kk.or.ii.eq.ll)goto 2000 if(jj.eq.kk.or.jj.eq.ll.or.kk.eq.ll)goto 2000 c call dangl(ii,jj,kk,ll,dan) <---- call edangl(ii,jj,kk,ll,dan,edan) <---- c WRITE(*,9)namat(ii),namat(jj),namat(kk),namat(ll),dan,edan, * strbuf(1:len_trim(strbuf)) WRITE(3,9)namat(ii),namat(jj),namat(kk),namat(ll),dan,edan, * strbuf(1:len_trim(strbuf)) goto 7 endif c C c------------------------------------------------------------------------------- c Equation of a plane c------------------------------------------------------------------------------- c if(ipar.eq.-1)then if(kk.lt.1.or.kk.gt.natoms)goto 2000 if(ii.eq.kk.or.jj.eq.kk)goto 2000 c call planeeq(ii,jj,kk,apl,dapl,bpl,dbpl,cpl,dcpl,dpl,ddpl) <---- c WRITE(*,19)'Equation of plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))), * apl,dapl, bpl,dbpl, cpl,dcpl, dpl,ddpl, * strbuf(1:len_trim(strbuf)) write(3,'(1x)') WRITE(3,19)'Equation of plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))), * apl,dapl, bpl,dbpl, cpl,dcpl, dpl,ddpl, * strbuf(1:len_trim(strbuf)) 19 format(1x,2a,', ',a,', ',a,' ] :'/ * 5x,'A = ',F20.10,' +-',F20.10,' * X'/ * 5x,'B = ',F20.10,' +-',F20.10,' * Y'/ * 5x,'C = ',F20.10,' +-',F20.10,' * Z'/ * 5x,'D = ',F20.10,' +-',F20.10,5x,a) c endif c C c------------------------------------------------------------------------------- c Angle between two planes c------------------------------------------------------------------------------- c if(ipar.eq.-2)then if(kk.lt.1.or.kk.gt.natoms)goto 2000 if(ll.lt.1.or.ll.gt.natoms)goto 2000 if(mm.lt.1.or.mm.gt.natoms)goto 2000 if(nn.lt.1.or.nn.gt.natoms)goto 2000 if(ii.eq.kk.or.jj.eq.kk)goto 2000 if(ll.eq.mm)goto 2000 if(ll.eq.nn.or.mm.eq.nn)goto 2000 c call plana(ii,jj,kk,ll,mm,nn,plangle) <---- call eplana(ii,jj,kk,ll,mm,nn,plangle,eplangle) <---- c WRITE(*,29)'Angle between plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))),' ]', * namat(ll)(1:len_trim(namat(ll))), * namat(mm)(1:len_trim(namat(mm))), * namat(nn)(1:len_trim(namat(nn))),' ]', * plangle,eplangle, * strbuf(1:len_trim(strbuf)) write(3,'(1x)') WRITE(3,29)'Angle between plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))),' ]', * namat(ll)(1:len_trim(namat(ll))), * namat(mm)(1:len_trim(namat(mm))), * namat(nn)(1:len_trim(namat(nn))),' ]', * plangle,eplangle, * strbuf(1:len_trim(strbuf)) 29 format(1x,2a,', ',a,', ',2a/10x,' and plane [ ', * a,', ',a,', ',2a/ * 25x,' = ',F10.5,' +-',F8.5,5x,a) c endif c C c------------------------------------------------------------------------------- c Angle between a plane and a line c------------------------------------------------------------------------------- c if(ipar.eq.-3)then if(kk.lt.1.or.kk.gt.natoms)goto 2000 if(ll.lt.1.or.ll.gt.natoms)goto 2000 if(mm.lt.1.or.mm.gt.natoms)goto 2000 if(ii.eq.kk.or.jj.eq.kk)goto 2000 if(ll.eq.mm)goto 2000 c call pllna(ii,jj,kk,ll,mm,pllinan) <---- call epllna(ii,jj,kk,ll,mm,pllinan,epllinan) <---- c WRITE(*,59)'Angle between plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))),' ]', * namat(ll)(1:len_trim(namat(ll))), * namat(mm)(1:len_trim(namat(mm))),' ]', * pllinan,epllinan, * strbuf(1:len_trim(strbuf)) write(3,'(1x)') WRITE(3,59)'Angle between plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))),' ]', * namat(ll)(1:len_trim(namat(ll))), * namat(mm)(1:len_trim(namat(mm))),' ]', * pllinan,epllinan, * strbuf(1:len_trim(strbuf)) 59 format(1x,2a,', ',a,', ',2a/10x,' and line [ ', * a,', ',2a/ * 25x,' = ',F10.5,' +-',F8.5,5x,a) c endif c C c------------------------------------------------------------------------------- c Angle between two lines c------------------------------------------------------------------------------- c if(ipar.eq.-4)then if(kk.lt.1.or.kk.gt.natoms)goto 2000 if(ll.lt.1.or.ll.gt.natoms)goto 2000 if(kk.eq.ll)goto 2000 c call lnlna(ii,jj,kk,ll,alinlin) <---- call elnlna(ii,jj,kk,ll,alinlin,ealinlin) <---- c WRITE(*,69)'Angle between line [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))),' ]', * namat(kk)(1:len_trim(namat(kk))), * namat(ll)(1:len_trim(namat(ll))),' ]', * alinlin,ealinlin, * strbuf(1:len_trim(strbuf)) write(3,'(1x)') WRITE(3,69)'Angle between line [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))),' ]', * namat(kk)(1:len_trim(namat(kk))), * namat(ll)(1:len_trim(namat(ll))),' ]', * alinlin,ealinlin, * strbuf(1:len_trim(strbuf)) 69 format(1x,2a,', ',2a/10x,' and line [ ', * a,', ',2a/ * 25x,' = ',F10.5,' +-',F8.5,5x,a) c endif c c c------------------------------------------------------------------------------- c Distance between a point and a plane c------------------------------------------------------------------------------- c if(ipar.eq.-5)then if(kk.lt.1.or.kk.gt.natoms)goto 2000 if(ll.lt.1.or.ll.gt.natoms)goto 2000 if(ii.eq.kk.or.jj.eq.kk)goto 2000 c call planeeq(ii,jj,kk,apl,dapl,bpl,dbpl,cpl,dcpl,dpl,ddpl) <---- ppdist=dabs( apl*cord(ll,1) * + bpl*cord(ll,2) * + cpl*cord(ll,3) * + dpl) / dsqrt(apl**2+bpl**2+cpl**2) c denom=dsqrt(apl**2+bpl**2+cpl**2) ddenom=dsqrt( (apl**2*dapl**2+bpl**2*dbpl**2+cpl**2*dcpl**2)/ * (apl**2+bpl**2+cpl**2) ) dnumer=dsqrt( cord(ll,1)**2*dapl**2 + apl**2*ecord(ll,1)**2 * + cord(ll,2)**2*dbpl**2 + bpl**2*ecord(ll,2)**2 * + cord(ll,3)**2*cbpl**2 + cpl**2*ecord(ll,3)**2 * + cdpl**2 ) dppdist=dsqrt(dnumer**2+ppdist**2*ddenom**2)/denom c WRITE(*,39)'Distance from plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))),' ]', * 'to point [ ', * namat(ll)(1:len_trim(namat(ll))), * ppdist,dppdist, * strbuf(1:len_trim(strbuf)) write(3,'(1x)') WRITE(3,39)'Distance from plane [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))), * namat(kk)(1:len_trim(namat(kk))),' ]', * 'to point [ ', * namat(ll)(1:len_trim(namat(ll))), * ppdist,dppdist, * strbuf(1:len_trim(strbuf)) 39 format(1x,2a,', ',a,', ',2a/ * 12x,2a,' ]'/ * 25x,' = ',F10.5,' +-',F8.5,5x,a) c endif c C c------------------------------------------------------------------------------- c Distance between a point and a line c------------------------------------------------------------------------------- c c parametric equation of a line through points P0 and P1 is C x = x0 + (x1-x0) Lambda C y = y0 + (y1-y0) Lambda C z = z0 + (z1-z0) Lambda c if(ipar.eq.-6)then aln=cord(jj,1)-cord(ii,1) bln=cord(jj,2)-cord(ii,2) cln=cord(jj,3)-cord(ii,3) daln=dsqrt(ecord(jj,1)**2+ecord(ii,1)**2) dbln=dsqrt(ecord(jj,2)**2+ecord(ii,2)**2) dcln=dsqrt(ecord(jj,3)**2+ecord(ii,3)**2) c ux=(cord(kk,1)-cord(ii,1)) uy=(cord(kk,2)-cord(ii,2)) uz=(cord(kk,3)-cord(ii,3)) dux=dsqrt(ecord(kk,1)**2+ecord(ii,1)**2) duy=dsqrt(ecord(kk,2)**2+ecord(ii,2)**2) duz=dsqrt(ecord(kk,3)**2+ecord(ii,3)**2) c aa=bln*uz-cln*uy bb=aln*uz-cln*ux cc=aln*uy-bln*ux daa=dsqrt( bln**2*duz**2 + uz**2*dbln**2 + * cln**2*duy**2 + uy**2*dcln**2 ) dbb=dsqrt( aln**2*duz**2 + uz**2*daln**2 + * cln**2*dux**2 + ux**2*dcln**2 ) dcc=dsqrt( aln**2*duy**2 + uy**2*daln**2 + * bln**2*dux**2 + ux**2*dbln**2 ) c top=dsqrt( aa**2+bb**2+cc**2) dtop=dsqrt( (aa**2*daa**2+bb**2*dbb**2+cc**2*dcc**2)/ * (aa**2+bb**2+cc**2) ) c denom=dsqrt(aln**2+bln**2+cln**2) ddenom=dsqrt( (aln**2*daln**2+bln**2*dbln**2+cln**2*dcln**2)/ * (aln**2+bln**2+cln**2) ) c pldist=top/denom dpldist=dsqrt(dtop**2+pldist**2*ddenom**2)/denom c WRITE(*,49)'Distance from line [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))),' ]', * 'to point [ ', * namat(ll)(1:len_trim(namat(kk))), * pldist,dpldist, * strbuf(1:len_trim(strbuf)) write(3,'(1x)') WRITE(3,49)'Distance from line [ ', * namat(ii)(1:len_trim(namat(ii))), * namat(jj)(1:len_trim(namat(jj))),' ]', * 'to point [ ', * namat(ll)(1:len_trim(namat(kk))), * pldist,dpldist, * strbuf(1:len_trim(strbuf)) 49 format(1x,2a,', ',a,a/ * 12x,2a,' ]'/ * 25x,' = ',F10.5,' +-',F8.5,5x,a) c endif c GOTO 7 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c 4 write(3,'(78(1h-))') close(3) write(*,2002)filout(1:len_trim(filout)) 2002 FORMAT(1X/' Output written to: 'a//) STOP c 2000 write(*,2001)errmes,nlinp,len_trim(line),line(1:len_trim(line)) 2001 Format(1x//' **** INPUT ERROR: ',a/20x, * 'when processing input line',i5/20x, * 'containing the following',i5,' characters:'/ * 1x,a//) close(3) c STOP END c C_____________________________________________________________________________ c subroutine bond(II,JJ,r) parameter (maxat=500) implicit real*8 (a-h,o-z) common /atoms/cord(maxat,3),ecord(maxat,3) a = cord(ii,1) b = cord(ii,2) c = cord(ii,3) a1 = cord(jj,1) b1 = cord(jj,2) c1 = cord(jj,3) DELA=A1-A DELB=B1-B DELC=C1-C r=dsqrt( DELA**2 + DELB**2 + DELC**2 ) C return end C C_____________________________________________________________________________ c subroutine eBOND(ii,jj,r,er) c c Numerical evaluation of the uncertainty in distance c implicit real*8 (a-h,o-z) parameter (maxat=500,cincr=0.001d0) real*8 cord(maxat,3),ecord(maxat,3),deriv(6,2),indx(6,2) common /atoms/cord,ecord c do 2 n=1,3 indx(n, 1)=ii indx(n+3,1)=jj indx(n, 2)=n indx(n+3,2)=n 2 continue c do 1 n=1,6 nn=indx(n,1) nnn=indx(n,2) cord(nn,nnn)=cord(nn,nnn)+cincr call bond(ii,jj,r1) <---- deriv(n,1)=(r1-r)/cincr deriv(n,2)=ecord(nn,nnn) cord(nn,nnn)=cord(nn,nnn)-cincr 1 continue c er=0.d0 do 3 n=1,6 er=er+deriv(n,1)**2*deriv(n,2)**2 3 continue er=dsqrt(er) c return end c C_____________________________________________________________________________ c subroutine ANGLE(ii,jj,kk,AN) c c The angle is that opposite r3 and is returned in degrees c parameter (maxat=500) implicit real*8 (a-h,o-z) common /atoms/cord(maxat,3),ecord(maxat,3) CONV=1.7453292519943D-2 2 pi/360 C call bond(ii,jj,r1) <---- Fortran2018 call bond(jj,kk,r2) <---- Fortran2018 call bond(ii,kk,r3) <---- Fortran2018 an=DACOS((r1**2+r2**2-r3**2)/(2.d0*r1*r2))/conv c RETURN END C C_____________________________________________________________________________ c subroutine eANGLE(ii,jj,kk,an,ean) c c Numerical evaluation of the uncertainty in an angle, which avoids c problems with analytical derivatives for small angles c implicit real*8 (a-h,o-z) parameter (maxat=500,cincr=0.001d0) real*8 cord(maxat,3),ecord(maxat,3),deriv(9,2),indx(9,2) common /atoms/cord,ecord c do 2 n=1,3 indx(n, 1)=ii indx(n+3,1)=jj indx(n+6,1)=kk indx(n, 2)=n indx(n+3,2)=n indx(n+6,2)=n 2 continue c do 1 n=1,9 nn=indx(n,1) nnn=indx(n,2) cord(nn,nnn)=cord(nn,nnn)+cincr call angle(ii,jj,kk,an1) <---- deriv(n,1)=(an1-an)/cincr deriv(n,2)=ecord(nn,nnn) cord(nn,nnn)=cord(nn,nnn)-cincr 1 continue c ean=0.d0 do 3 n=1,9 ean=ean+deriv(n,1)**2*deriv(n,2)**2 3 continue ean=dsqrt(ean) c return end c C_____________________________________________________________________________ c subroutine DANGL(ii,jj,kk,ll,dan) C C Dihedral angle (ii.jj.kk.ll), returned in degrees C parameter (maxat=500) implicit real*8 (a-h,o-z) real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord CONV=1.7453292519943D-2 2 pi/360 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 call angle(ii,jj,kk,ANGL1) <----- Fortran2018 ANGL1=ANGL1-90.d0 call angle(jj,kk,ll,ANGL2) <----- fortran2018 ANGL2=ANGL2-90.d0 c call bond(ii,jj,R12) <----- Fortran2018 R12S=R12*DSIN(ANGL1*CONV) c call bond(jj,kk,R23) <----- Fortran2018 call bond(kk,ll,R34) <----- Fortran2018 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) C IF(DABS(DAN).GT.1.D0)DAN=DSIGN(1.D0,DAN) DAN=DACOS(DAN)/CONV C C Sign check (by Marek Jaworski, 18.11.2011) C aa1=a1-b1 aa2=a2-b2 aa3=a3-b3 cc1=c1-b1 cc2=c2-b2 cc3=c3-b3 dd1=d1-b1 dd2=d2-b2 dd3=d3-b3 det=cc1*(aa2*dd3-aa3*dd2)+cc2*(aa3*dd1-aa1*dd3) det=det+cc3*(aa1*dd2-aa2*dd1) if(det.lt.0d0)dan=-dan C RETURN END c C_____________________________________________________________________________ c subroutine eDANGL(ii,jj,kk,ll,dan,edan) c c Numerical evaluation of the uncertainty in the dihedral angle, which avoids c problems with analytical derivatives for small angles c implicit real*8 (a-h,o-z) parameter (maxat=500,cincr=0.001d0) real*8 cord(maxat,3),ecord(maxat,3),deriv(12,2),indx(12,2) common /atoms/cord,ecord c do 2 n=1,3 indx(n, 1)=ii indx(n+3,1)=jj indx(n+6,1)=kk indx(n+9,1)=ll indx(n, 2)=n indx(n+3,2)=n indx(n+6,2)=n indx(n+9,2)=n 2 continue c do 1 n=1,12 nn=indx(n,1) nnn=indx(n,2) cord(nn,nnn)=cord(nn,nnn)+cincr call dangl(ii,jj,kk,ll,dan1) <---- deriv(n,1)=(dan1-dan)/cincr deriv(n,2)=ecord(nn,nnn) cord(nn,nnn)=cord(nn,nnn)-cincr 1 continue c edan=0.d0 do 3 n=1,12 edan=edan+deriv(n,1)**2*deriv(n,2)**2 3 continue edan=dsqrt(edan) c return end c C_____________________________________________________________________________ c subroutine planeeq(ii,jj,kk,apl,dapl,bpl,dbpl,cpl,dcpl,dpl,ddpl) c c Equation of a plane containing points ii,jj,kk: c APL*x + BPL*y + CPL*z + DPL = 0 c implicit real*8 (a-h,o-z) parameter (maxat=500,cincr=0.001d0) real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord c c...determine vectors v=ii-jj and u=ii-kk c vx=cord(ii,1)-cord(jj,1) dvx=dsqrt(ecord(ii,1)**2+ecord(jj,1)**2) vy=cord(ii,2)-cord(jj,2) dvy=dsqrt(ecord(ii,2)**2+ecord(jj,2)**2) vz=cord(ii,3)-cord(jj,3) dvz=dsqrt(ecord(ii,3)**2+ecord(jj,3)**2) ux=cord(ii,1)-cord(kk,1) dux=dsqrt(ecord(ii,1)**2+ecord(kk,1)**2) uy=cord(ii,2)-cord(kk,2) duy=dsqrt(ecord(ii,2)**2+ecord(kk,2)**2) uz=cord(ii,3)-cord(kk,3) duz=dsqrt(ecord(ii,3)**2+ecord(kk,3)**2) c c...normal to the plane given by the cross product c v x u = det | v1 v2 v3 | c | u1 u2 u3 | c apl= (vy*uz-uy*vz) dapl=dsqrt( uz**2*dvy**2 + vy**2*duz**2 * + vz**2*duy**2 + uy**2*dvz**2) bpl=-(vx*uz-ux*vz) dbpl=dsqrt( uz**2*dvx**2 + vx**2*duz**2 * + vz**2*dux**2 + ux**2*dvz**2) cpl= (vx*uy-ux*vy) dcpl=dsqrt( uy**2*dvx**2 + vx**2*duy**2 * + vy**2*dux**2 + ux**2*dvy**2) dpl=-(apl*cord(ii,1)+bpl*cord(ii,2)+cpl*cord(ii,3)) ddpl=dsqrt( cord(ii,1)**2*dapl**2 + apl**2*ecord(ii,2)**2 * + cord(ii,2)**2*dbpl**2 + bpl**2*ecord(ii,2)**2 * + cord(ii,3)**2*dcpl**2 + cpl**2*ecord(ii,3)**2 ) c return end c C_____________________________________________________________________________ c subroutine plana(ii,jj,kk,ll,mm,nn,plangle) C C Angle between planes (ii.jj.kk) and (ll,mm,nn) returned in degrees C parameter (maxat=500) implicit real*8 (a-h,o-z) real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord CONV=1.7453292519943D-2 2 pi/360 c call planeeq(ii,jj,kk,apl,dapl,bpl,dbpl,cpl,dcpl,dpl,ddpl) <---- call planeeq(ll,mm,nn,a1pl,da1pl,b1pl,db1pl, * c1pl,dc1pl,d1pl,dd1pl) <---- c top=dabs(apl*a1pl+bpl*b1pl+cpl*c1pl) c denom1=dsqrt(apl**2+bpl**2+cpl**2) denom2=dsqrt(a1pl**2+b1pl**2+c1pl**2) denom=denom1*denom2 vcos=top/denom plangle=dacos( vcos )/conv c return end c C_____________________________________________________________________________ c subroutine eplana(ii,jj,kk,ll,mm,nna,plangle,eplangle) c c Numerical evaluation of the uncertainty in the interplane angle, which avoids c problems with analytical derivatives for small angles c implicit real*8 (a-h,o-z) parameter (maxat=500,cincr=0.001d0) real*8 cord(maxat,3),ecord(maxat,3),deriv(18,2),indx(18,2) common /atoms/cord,ecord c do 2 n=1,3 indx(n, 1)=ii indx(n+ 3,1)=jj indx(n+ 6,1)=kk indx(n+ 9,1)=ll indx(n+12,1)=mm indx(n+15,1)=nna indx(n, 2)=n indx(n+ 3,2)=n indx(n+ 6,2)=n indx(n+ 9,2)=n indx(n+12,2)=n indx(n+15,2)=n 2 continue c do 1 n=1,18 nn=indx(n,1) nnn=indx(n,2) cord(nn,nnn)=cord(nn,nnn)+cincr call plana(ii,jj,kk,ll,mm,nna,plangle1) <---- deriv(n,1)=(plangle1-plangle)/cincr deriv(n,2)=ecord(nn,nnn) cord(nn,nnn)=cord(nn,nnn)-cincr 1 continue c eplangle=0.d0 do 3 n=1,18 eplangle=eplangle+deriv(n,1)**2*deriv(n,2)**2 3 continue eplangle=dsqrt(eplangle) c return end c C_____________________________________________________________________________ c subroutine pllna(ii,jj,kk,ll,mm,pllinan) c C Angle between plane (ii.jj.kk) and line (ll,mm) returned in degrees C parameter (maxat=500) implicit real*8 (a-h,o-z) real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord CONV=1.7453292519943D-2 2 pi/360 c call planeeq(ii,jj,kk,apl,dapl,bpl,dbpl,cpl,dcpl,dpl,ddpl) <---- c ux=cord(ll,1)-cord(mm,1) uy=cord(ll,2)-cord(mm,2) uz=cord(ll,3)-cord(mm,3) c top=dabs( ux*apl+uy*bpl+uz*cpl ) denom=dsqrt( apl**2+bpl**2+cpl**2)* * dsqrt( ux**2+uy**2+uz**2) c pllinan=top/denom pllinan=dasin(pllinan)/conv c return end c C_____________________________________________________________________________ c subroutine epllna(ii,jj,kk,ll,mm,pllinan,epllinan) c c Numerical evaluation of the uncertainty in the plane-line angle, which avoids c problems with analytical derivatives for small angles c implicit real*8 (a-h,o-z) parameter (maxat=500,cincr=0.001d0) real*8 cord(maxat,3),ecord(maxat,3),deriv(15,2),indx(15,2) common /atoms/cord,ecord c do 2 n=1,3 indx(n, 1)=ii indx(n+ 3,1)=jj indx(n+ 6,1)=kk indx(n+ 9,1)=ll indx(n+12,1)=mm indx(n, 2)=n indx(n+ 3,2)=n indx(n+ 6,2)=n indx(n+ 9,2)=n indx(n+12,2)=n 2 continue c do 1 n=1,15 nn=indx(n,1) nnn=indx(n,2) cord(nn,nnn)=cord(nn,nnn)+cincr call pllna(ii,jj,kk,ll,mm,pllinan1) <---- deriv(n,1)=(pllinan1-pllinan)/cincr deriv(n,2)=ecord(nn,nnn) cord(nn,nnn)=cord(nn,nnn)-cincr 1 continue c epllinan=0.d0 do 3 n=1,15 Fortran2018 epllinan=epllinan+deriv(n,1)**2*deriv(n,2)**2 3 continue epllinan=dsqrt(epllinan) c return end c C_____________________________________________________________________________ c subroutine lnlna(ii,jj,kk,ll,alinlin) c C Angle between line (ii.jj) and line (kk,ll) returned in degrees C parameter (maxat=500) implicit real*8 (a-h,o-z) real*8 cord(maxat,3),ecord(maxat,3) common /atoms/cord,ecord CONV=1.7453292519943D-2 2 pi/360 c ux=cord(ii,1)-cord(jj,1) uy=cord(ii,2)-cord(jj,2) uz=cord(ii,3)-cord(jj,3) vx=cord(kk,1)-cord(ll,1) vy=cord(kk,2)-cord(ll,2) vz=cord(kk,3)-cord(ll,3) c top=dabs( ux*vx+uy*vy+uz*vz ) denom=dsqrt( ux**2+uy**2+uz**2)* * dsqrt( vx**2+vy**2+vz**2) c alinlin=top/denom alinlin=dacos(alinlin)/conv c return end c C_____________________________________________________________________________ c subroutine elnlna(ii,jj,kk,ll,alinlin,ealinlin) c c Numerical evaluation of the uncertainty in the line-line angle, which avoids c problems with analytical derivatives for small angles c implicit real*8 (a-h,o-z) parameter (maxat=500,cincr=0.001d0) real*8 cord(maxat,3),ecord(maxat,3),deriv(12,2),indx(12,2) common /atoms/cord,ecord c do 2 n=1,3 indx(n, 1)=ii indx(n+ 3,1)=jj indx(n+ 6,1)=kk indx(n+ 9,1)=ll indx(n, 2)=n indx(n+ 3,2)=n indx(n+ 6,2)=n indx(n+ 9,2)=n 2 continue c do 1 n=1,12 nn=indx(n,1) nnn=indx(n,2) cord(nn,nnn)=cord(nn,nnn)+cincr call lnlna(ii,jj,kk,ll,alinlin1) <---- deriv(n,1)=(alinlin1-alinlin)/cincr deriv(n,2)=ecord(nn,nnn) cord(nn,nnn)=cord(nn,nnn)-cincr 1 continue c ealinlin=0.d0 do 3 n=1,12 Fortran2018 ealinlin=ealinlin+deriv(n,1)**2*deriv(n,2)**2 3 continue ealinlin=dsqrt(ealinlin) c return end c C_____________________________________________________________________________ c subroutine test(iodev,inperr) character line*79 c c...test whether the next line in currently open file on IODEV is a comment c line (INPERR=1), otherwise INPERR=0 c read(iodev,'(a)',err=3,end=3)line if(line(1:1).eq.'!'.or.line(1:1).eq.'#'.or.line(1:1).eq.'%')then inperr=1 if(line(1:1).eq.'!'.and.line(1:6).ne.'!NNNNN')then write(*,'(1x,a)')line(1:len_trim(line)) write(3,'(1x,a)')line(1:len_trim(line)) endif else backspace(iodev) inperr=0 endif return c 3 inperr=2 write(*,*)' --- end of input file ' return end c C_____________________________________________________________________________ C_____________________________________________________________________________