c23456789012345678901234567890123456789012345678901234567890123456789012
program convpdbfi
character*3 resnam
character*3 restyp(800)
character*8 fil(160)
character*4 prna
character*4 atnam
character*6 dummy
integer atnum,resnum,resold,maxres,icount(22)
integer nres(800),natoms(800),nala(22,800),iala(22,800),iss2(200)
integer dist(22,25),resno(160),iss1(200)
integer icounf(22,22),totfi(22,22),ibond(22)
integer icounk(22),iir(150,800),doublet(22,13,13)
integer dt(22,22,18),isumxr(22),isumrx(22),doublet2(22,13,13)
integer doublet3(22,13,13),dall(13,13),resno1(800)
real x(800,35),y(800,35),z(800,35),totan(22),avgan(22)
real mass(800,35),xx(800),yy(800),zz(800)
real avgbo(22),ax(6,80,2),ay(6,80,2),az(6,80,2)
real elx(800),ely(800),elz(800),elm(800),elmsq(800)
real crox(800),croy(800),croz(800),crom(800),cromsq(800)
real croxp(800),croyp(800),crozp(800),cromp(800),hyd(21)
real theta(150,800),phi(150,800)
real bond(22,4700),thet1(22),thet2(22)
real sidex(800),sidey(800),sidez(800)
real cyss(150,800),blength(150,800)
real thetpr(150,800),phip(150,800)
c This program calculates virtual bond model coordinates. First
c virtual bond lengths are evaluated. Then bond angles will be
c determined for each residue. Finally torsional angles between
c virtual bonds will be evaluated. We need a single point per
c residue, which will be found from the mass center of the points
c already considered in evaluating pair radial distribution fns.
c These points will be saved as xx, yy and zz, for each residue in
c a given protein. Therefore the corresponding array sizes will
c be 800. Likewise the alpha carbons will be used for specifying
c backbone locations
c nala(kk,iiii) designates the number of residues of type kk in
c the protein iiii. natoms(kk) designates the total number of
c atoms of residue kk, in a given protein, whose coordinates are
c listed in PDB. iala(kk,i) denotes the serial index of the ith
c residue of type kk along a given protein.
fil(1)='1STN.pdb'
do 59331 i=1,22
icount(i)=0
ibond(i)=0
do 59332 j=1,22
59332 icounf(i,j)=0
59331 continue
c nfile2>1 if more than 1 file is to be read
nfile1=1
nfile2=1
c_____________________________________________________________________
do 4293 iiii=nfile1,nfile2
do 40893 j=1,800
xx(j)=0.
yy(j)=0.
40893 zz(j)=0.
nsb=0
conv=180./3.141592654
write(6,*)'********',iiii
open(10,file=fil(iiii))
do 3 i=1,600
do 3 j=1,35
x(i,j)=0.
y(i,j)=0.
3 z(i,j)=0.
i=0
resold=1
5 read(10,50)dummy
50 Format(A6)
40 format (A4)
55 FORMAT (A6,1X,I4,1X,A4,1X,A3,1X,A1,1X,I3,5X,3F8.3,5X,F6.2)
54 format (A6,1X,I3,1X,A1,2X,I3,2X,13(A3,1X),1X,A4)
56 format (A6,1X,I3,1X,A3,1X,A1,2X,I3,4X,A3,1X,A1,2X,I3)
if (dummy.eq.'SEQRES') then
backspace (10)
read(10,54)dummy,nrow,s,maxres,(restyp(13*(nrow-1)+j),j=1,13),prna
write(70,54)dummy,nrow,s,maxres,(restyp(13*(nrow-1)+j),j=1,13),
+ prna
endif
maxres=maxresold
if(dummy.ne.'ATOM ') goto 5
backspace (10)
12 i=i+1
1231 read(10,50)dummy
if(dummy.eq.'HETATM') goto 33
if(dummy.eq.'SIGATM') goto 1231
backspace (10)
read(10,55)dummy,atnum,atnam,resnam,s,resnum,x(resnum,i),
+ y(resnum,i),z(resnum,i),bf
write(70,55)dummy,atnum,atnam,resnam,s,resnum,
+ x(resnum,i),y(resnum,i),z(resnum,i),bf
if(atnam.eq.' CA ') then
write(92,*)resnum,bf
endif
if(atnum.eq.1) resno1(iiii)=resnum
if(resnum.ne.resold)then
natoms(resold)=i-1
x(resnum,1)=x(resnum,i)
y(resnum,1)=y(resnum,i)
z(resnum,1)=z(resnum,i)
restyp(resnum)=resnam
i=1
resold=resnum
endif
if(dummy.eq.'TER ') goto 33
c if(atnum.eq.757) goto 33
go to 12
c*************************************************************
33 maxat=atnum
c write(91,*)'resnum',resnum
natoms(resnum)=i-1
resno(iiii)=resnum
do 29811 i=1,resnum
blength(iiii,i)=0.
theta(iiii,i)=0.
thetpr(iiii,i)=0.
phip(iiii,i)=0.
29811 phi(iiii,i)=0.
kala=1
kgly=1
kval=1
kleu=1
kile=1
kphe=1
karg=1
klys=1
kglu=1
kasp=1
kasn=1
kgln=1
kpro=1
khis=1
ktyr=1
ktrp=1
kthr=1
kser=1
kcys=1
kmet=1
kcyss=1
do 19001 i=resno1(iiii),resnum
write(91,*)i,restyp(i)
c atomic mass of different atoms N, O, C..
mass(i,1)=15.
mass(i,2)=13.
mass(i,3)=12.
mass(i,4)=16.
do 12022 j=5,14
12022 mass(i,j)=12.
do 12029 j=15,25
12029 mass(i,j)=1.
if (restyp(i).eq.'ALA') THEN
iir(iiii,i)=2
if (natoms(i).ge.5) then
write(91,*)'kkkkkkk',i
iala(2,kala)=i
kala=kala+1
xx(i)=x(i,5)
yy(i)=y(i,5)
zz(i)=z(i,5)
endif
endif
if (restyp(i).eq.'VAL') THEN
iir(iiii,i)=3
if (natoms(i).ge.7) then
write(91,*)'kkkkkkk',i
iala(3,kval)=i
kval=kval+1
xx(i)=(x(i,6)+x(i,7))/2.0
yy(i)=(y(i,6)+y(i,7))/2.0
zz(i)=(z(i,6)+z(i,7))/2.0
endif
endif
if (restyp(i).eq.'LEU') THEN
iir(iiii,i)=5
if (natoms(i).ge.8) then
write(91,*)'kkkkkkk',i
iala(5,kleu)=i
kleu=kleu+1
xx(i)=(x(i,8)+x(i,7))/2.0
yy(i)=(y(i,8)+y(i,7))/2.0
zz(i)=(z(i,8)+z(i,7))/2.0
endif
endif
if (restyp(i).eq.'ILE') THEN
iir(iiii,i)=4
if (natoms(i).ge.8) then
write(91,*)'kkkkkkk',i
iala(4,kile)=i
kile=kile+1
xx(i)=x(i,8)
yy(i)=y(i,8)
zz(i)=z(i,8)
endif
endif
if (restyp(i).eq.'PHE') THEN
iir(iiii,i)=16
mass(i,2)=17.
if (natoms(i).ge.11) then
write(91,*)'kkkkkkk',i
iala(16,kphe)=i
kphe=kphe+1
do 48331 iji=6,11
xx(i)=xx(i)+x(i,iji)
yy(i)=yy(i)+y(i,iji)
48331 zz(i)=zz(i)+z(i,iji)
xx(i)=xx(i)/6.
yy(i)=yy(i)/6.
zz(i)=zz(i)/6.
endif
endif
if (restyp(i).eq.'CYS') THEN
do 38999 ij=1,nsb
if(i.eq.iss1(ij).or.i.eq.iss2(ij))then
iala(21,kcyss)=i
cyss(iiii,i)=1.
kcyss=kcyss+1
endif
38999 continue
iala(14,kcys)=i
kcys=kcys+1
mass(i,6)=33.
iir(iiii,i)=14
if (natoms(i).ge.6) then
write(91,*)'kkkkkkk',i
yy(i)=y(i,6)
zz(i)=z(i,6)
xx(i)=x(i,6)
endif
endif
if (restyp(i).eq.'MET') THEN
iir(iiii,i)=15
mass(i,7)=32.
if (natoms(i).ge.7) then
write(91,*)'kkkkkkk',i
iala(15,kmet)=i
kmet=kmet+1
xx(i)=x(i,7)
yy(i)=y(i,7)
zz(i)=z(i,7)
endif
endif
if (restyp(i).eq.'HIS') then
iir(iiii,i)=19
mass(i,7)=14.
mass(i,10)=14.
if (natoms(i).ge.10) then
write(91,*)'kkkkkkk',i
iala(19,khis)=i
khis=khis+1
do 48332 iji=6,10
xx(i)=xx(i)+x(i,iji)
yy(i)=yy(i)+y(i,iji)
48332 zz(i)=zz(i)+z(i,iji)
xx(i)=xx(i)/5.
yy(i)=yy(i)/5.
zz(i)=zz(i)/5.
endif
endif
if (restyp(i).eq.'TRP') then
iir(iiii,i)=18
mass(i,9)=14.
if (natoms(i).ge.14) then
write(91,*)'kkkkkkk',i
iala(18,ktrp)=i
ktrp=ktrp+1
do 48333 iji=7,14
xx(i)=xx(i)+x(i,iji)
yy(i)=yy(i)+y(i,iji)
48333 zz(i)=zz(i)+z(i,iji)
xx(i)=xx(i)/8.
yy(i)=yy(i)/8.
zz(i)=zz(i)/8.
endif
endif
if (restyp(i).eq.'ASP') then
iir(iiii,i)=8
mass(i,7)=16.
mass(i,8)=16.
if (natoms(i).ge.8) then
write(91,*)'kkkkkkk',i
iala(8,kasp)=i
kasp=kasp+1
xx(i)=(x(i,8)+x(i,7))/2.0
yy(i)=(y(i,8)+y(i,7))/2.0
zz(i)=(z(i,8)+z(i,7))/2.0
endif
endif
if (restyp(i).eq.'GLU') then
iir(iiii,i)=10
mass(i,9)=16.
mass(i,8)=16.
if (natoms(i).ge.9) then
write(91,*)'kkkkkkk',i
iala(10,kglu)=i
kglu=kglu+1
xx(i)=(x(i,8)+x(i,9))/2.0
yy(i)=(y(i,8)+y(i,9))/2.0
zz(i)=(z(i,8)+z(i,9))/2.0
endif
endif
if (restyp(i).eq.'ASN') then
iir(iiii,i)=9
mass(i,7)=16.
mass(i,8)=14.
if (natoms(i).ge.8) then
write(91,*)'kkkkkkk',i
iala(9,kasn)=i
kasn=kasn+1
xx(i)=(x(i,8)+x(i,7))/2.0
yy(i)=(y(i,8)+y(i,7))/2.0
zz(i)=(z(i,8)+z(i,7))/2.0
endif
endif
if (restyp(i).eq.'GLN') then
iir(iiii,i)=11
mass(i,8)=16.
mass(i,9)=14.
if (natoms(i).ge.9) then
write(91,*)'kkkkkkk',i
iala(11,kgln)=i
kgln=kgln+1
xx(i)=(x(i,8)+x(i,9))/2.0
yy(i)=(y(i,8)+y(i,9))/2.0
zz(i)=(z(i,8)+z(i,9))/2.0
endif
endif
if (restyp(i).eq.'SER') then
iir(iiii,i)=6
mass(i,6)=16.
if (natoms(i).ge.6) then
write(91,*)'kkkkkkk',i
iala(6,kser)=i
kser=kser+1
xx(i)=x(i,6)
yy(i)=y(i,6)
zz(i)=z(i,6)
endif
endif
if (restyp(i).eq.'THR') then
iir(iiii,i)=7
mass(i,6)=16.
if (natoms(i).ge.6) then
write(91,*)'kkkkkkk',i
iala(7,kthr)=i
kthr=kthr+1
xx(i)=x(i,6)
yy(i)=y(i,6)
zz(i)=z(i,6)
endif
endif
if (restyp(i).eq.'LYS') then
iir(iiii,i)=12
mass(i,9)=15.
if (natoms(i).ge.9) then
write(91,*)'kkkkkkk',i
iala(12,klys)=i
klys=klys+1
xx(i)=x(i,9)
yy(i)=y(i,9)
zz(i)=z(i,9)
endif
endif
if (restyp(i).eq.'ARG') then
iir(iiii,i)=13
mass(i,8)=15.
mass(i,10)=16.
mass(i,11)=16.
if (natoms(i).ge.11) then
write(91,*)'kkkkkkk',i
iala(13,karg)=i
karg=karg+1
xx(i)=(x(i,8)+x(i,10)+x(i,11))/3.0
yy(i)=(y(i,8)+y(i,10)+y(i,11))/3.0
zz(i)=(z(i,8)+z(i,10)+z(i,11))/3.0
endif
endif
if (restyp(i).eq.'TYR') then
iir(iiii,i)=17
mass(i,12)=17.
if (natoms(i).ge.12) then
write(91,*)'kkkkkkk',i
iala(17,ktyr)=i
ktyr=ktyr+1
do 48334 iji=6,12
xx(i)=xx(i)+x(i,iji)
yy(i)=yy(i)+y(i,iji)
48334 zz(i)=zz(i)+z(i,iji)
xx(i)=xx(i)/7.
yy(i)=yy(i)/7.
zz(i)=zz(i)/7.
endif
endif
if (restyp(i).eq.'GLY') then
iir(iiii,i)=1
if (natoms(i).ge.1) then
write(91,*)'kkkkkkk',i
iala(1,kgly)=i
kgly=kgly+1
xx(i)=x(i,2)
yy(i)=y(i,2)
zz(i)=z(i,2)
endif
endif
if (restyp(i).eq.'PRO') then
iir(iiii,i)=20
if (natoms(i).ge.7) then
write(91,*)'kkkkkkk',i
iala(20,kpro)=i
kpro=kpro+1
xx(i)=(x(i,5)+x(i,6)+x(i,7))/3.0
yy(i)=(y(i,5)+y(i,6)+y(i,7))/3.0
zz(i)=(z(i,5)+z(i,6)+z(i,7))/3.0
endif
endif
19001 continue
sum=0.
sum=kala+kgly+kval+kleu+kile+kphe+karg+klys+kglu+kasp+
+ kasn+kgln+kpro+khis+ktyr+ktrp+kthr+kser+kcys+kmet
c write(91,*)'sum',sum
kala=kala-1
kgly=kgly-1
kval=kval-1
kleu=kleu-1
kile=kile-1
kphe=kphe-1
karg=karg-1
klys=klys-1
kglu=kglu-1
kasp=kasp-1
kasn=kasn-1
kgln=kgln-1
kpro=kpro-1
khis=khis-1
ktyr=ktyr-1
ktrp=ktrp-1
kthr=kthr-1
kser=kser-1
kcys=kcys-1
kmet=kmet-1
kcyss=kcyss-1
nala(2,iiii)=kala
nala(1,iiii)=kgly
nala(3,iiii)=kval
nala(4,iiii)=kile
nala(5,iiii)=kleu
nala(6,iiii)=kser
nala(7,iiii)=kthr
nala(8,iiii)=kasp
nala(9,iiii)=kasn
nala(10,iiii)=kglu
nala(11,iiii)=kgln
nala(12,iiii)=klys
nala(13,iiii)=karg
nala(14,iiii)=kcys
nala(15,iiii)=kmet
nala(16,iiii)=kphe
nala(17,iiii)=ktyr
nala(18,iiii)=ktrp
nala(19,iiii)=khis
nala(20,iiii)=kpro
nala(21,iiii)=kcyss
sum=0
do 89 k=1,20
sum=sum+nala(k,1)
89 continue
write(*,*)'number of sum',sum,'unrecorded',maxres-sum
c__________________________________________________________________
c Residues with unknown coords are excluded from calculations
c__________________________________________________________________
c CALCULATION OF BACKBONE BOND LENGTHS AND ANGLES
DO 50572 K=resno1(iiii)+1,resnum
CROX(K) = 0.
CROY(K) = 0.
CROZ(K) = 0.
ELX(K) = 0.
ELY(K) = 0.
50572 ELZ(K) = 0.
DO 50571 K=resno1(iiii)+1,resnum
if (x(K,2).eq.0.or.x(K-1,2).eq.0)then
write(91,*)'kkkkkkkkk',k,iir(iiii,k)
goto 50571
endif
ELX(K) = x(K,2) - x(K-1,2)
ELY(K) = y(K,2) - y(K-1,2)
ELZ(K) = z(K,2) - z(K-1,2)
ELMSQ(K) = ELX(K)*ELX(K) + ELY(K)*ELY(K) + ELZ(K)*ELZ(K)
write(91,*)sqrt(ELMSQ(K)),k
50571 ELM(K) = SQRT (ELMSQ(K))
c write(6,*)iiii,elm(10),elm(20),elm(30),elm(40)
if (x(3,2).eq.0.or.x(2,2).eq.0) goto 50576
ELX(1) = x(3,2) - x(2,2)
ELY(1) = y(3,2) - y(2,2)
ELZ(1) = z(3,2) - z(2,2)
ELMSQ(1) = ELX(1)*ELX(1) + ELY(1)*ELY(1) + ELZ(1)*ELZ(1)
ELM(1) = SQRT (ELMSQ(1))
50576 k=resnum-1
kf=resnum+1
if (x(k,2).eq.0.or.x(k-1,2).eq.0) goto 50577
ELX(kf) = x(k,2) - x(k-1,2)
ELY(kf) = y(k,2) - y(k-1,2)
ELZ(kf) = z(k,2) - z(k-1,2)
ELMSQ(kf) = ELX(kf)*ELX(kf) + ELY(kf)*ELY(kf) + ELZ(kf)*ELZ(kf)
ELM(kf) = SQRT (ELMSQ(kf))
50577 continue
c__________________________________________________________________
c write(2,*)'CALCULATION OF VIRTUAL BOND ANGLES'
c WRITE(2,*)' '
c here theta(iiii,mi) means the bond angle of the mi th residue
c in the iiii th protein.phi(iiii,k)
c has similar meaning as theta(iiii,mi).
do 10000 jres1=1,20
DO 40571 i=1,nala(jres1,iiii)
mi=iala(jres1,i)
if (mi.lt.2) goto 40571
if (mi.gt.(resnum-1)) goto 40571
if (elx(mi).eq.0.or.elx(mi+1).eq.0) then
write(6,*)'one zero',iiii,mi,jres1
goto 40571
endif
icount(jres1)=icount(jres1)+1
kc=icount(jres1)
dot=elx(mi+1)*elx(mi)+ely(mi+1)*ely(mi)
dot = -(dot + elz(mi+1)*elz(mi))/elm(mi)/elm(mi+1)
theta(iiii,mi)=acos(dot)*conv
write(91,*)'th',theta(iiii,mi),mi
40571 continue
totan(jres1)=icount(jres1)
10000 continue
c__________________________________________________________________
C **CROSS PRODUCTS OF BOND VECTORS I AND (I+1) AND MAGNITUDES**
c here we include the hypothetical first and last bonds as well.
DO 801 K = resno1(iiii),resnum
if (elx(k).ne.0.and.elx(k+1).ne.0) then
croX(K) = ELY(K)*ELZ(K+1)-ELZ(K)*ELY(K+1)
croY(K) = -ELX(K)*ELZ(K+1)+ELZ(K)*ELX(K+1)
croZ(K) = ELX(K)*ELY(K+1)-ELY(K)*ELX(K+1)
croM(K) = crox(K)*crox(k)+croy(K)*croy(k)+croz(K)*croz(k)
crom(k)=crom(k)**0.5
c write(6,*)elx(k),ely(k),elz(k),elx(k+1),ely(k+1),elz(k+1)
c write(6,*)'cros pr',iiii,k,crox(k),croy(k),croz(k),crom(k)
endif
801 CONTINUE
c write(2,*)'CALCULATION OF TORSIONAL ANGLES'
DO 804 K = resno1(iiii),resnum
if (crox(k).ne.0.and.crox(k-1).ne.0) then
ara=CROX(K)*CROX(K-1)+CROY(K)*CRoy(k-1)+CROZ(K)*CROZ(K-1)
EKFY = -1./CROM(K)/CROM(K-1)
cosfi = ara*EKFY
if (cosfi.gt.1.) cosfi=1.0
if (cosfi.lt.-1.) cosfi=-1.0
FI = ACOS(cosfi)
gpm = crox(k-1)*elx(k+1)+croy(k-1)*ely(k+1)+croz(k-1)*elz(k+1)
if (gpm.gt.0.0) fi=-fi
faydeg = fi*conv
c write(6,*)'angle',iiii,k,faydeg,cosfi,ara
i1=iir(iiii,k-1)
i2=iir(iiii,k)
icounf(i1,i2)=icounf(i1,i2)+1
kc=icounf(i1,i2)
phi(iiii,k)=faydeg
c write(6,*)' '
write(91,*)'fay',phi(iiii,k),k
endif
phi(iiii,2)=0.
804 CONTINUE
c__________________________________________________________________
c Now we know each bond bending, torsion, location of side groups
c in space. Locations wrt backbone will be described by three
c variables: these are one length and two orientations as follows:
c (1)length of side group (i.e.distance between alpha-C
c and center xx,yy,zz of side group), (2) The angle the side group
c i makes with the plane of the alpha carbons i-1, i and i+1,
c and (3) The angle the side group makes with the normal plane to
c bisecting the bond angle i.
do 10032 k=resno1(iiii),resnum
if (restyp(k).eq.'GLY') then
sidex(k)=0.
sidey(k)=0.
sidez(k)=0.
goto 10032
endif
if (xx(k).ne.0.0) then
x1=xx(k)-x(k,2)
y1=yy(k)-y(k,2)
z1=zz(k)-z(k,2)
ara = x1*x1 + y1*y1 + z1*z1
ara=SQRT(ara)
blength(iiii,k)=ara
jres=iir(iiii,k)
ibond(jres)=ibond(jres)+1
mco=ibond(jres)
bond(jres,mco)=ara
sidex(k)=x1/bond(jres,mco)
sidey(k)=y1/bond(jres,mco)
sidez(k)=z1/bond(jres,mco)
c write(6,*)iiii,k,jres,ibond(jres),bond(jres,mco)
endif
10032 continue
do 10001 jres1=1,20
DO 42571 i=1,nala(jres1,iiii)
mi=iala(jres1,i)
if (restyp(mi).eq.'GLY') THEN
thetpr(iiii,mi)=0.
phip(iiii,mi)=0.
goto 42571
endif
if (elx(mi).eq.0.or.sidex(mi).eq.0) then
write(6,*)'one side group zero',iiii,mi,jres1
goto 42571
endif
dot=sidex(mi)*elx(mi)+sidey(mi)*ely(mi)
dot = -(dot + sidez(mi)*elz(mi))/elm(mi)
thetpr(iiii,mi)=acos(dot)*conv
42571 continue
10001 continue
if (elx(1).eq.0.or.sidex(1).eq.0) then
write(6,*)'one terminal side group zero'
goto 42574
endif
dot=-sidex(1)*elx(2)-sidey(1)*ely(2)
dot = -(dot - sidez(1)*elz(2))/elm(2)
thetpr(iiii,1)=acos(dot)*conv
42574 continue
c cross product of bond i and side bond appended to residue i
DO 80117 K = resno1(iiii)+1,resnum
if (elx(k).ne.0.and.sidex(k).ne.0) then
croxp(K) = ELY(K)*sideZ(K)-ELZ(K)*sideY(K)
c write(6,*)'now',croxp(k),ely(k),elz(k),sidez(k),sidey(k)
croyp(K) = -ELX(K)*sideZ(K)+ELZ(K)*sideX(K)
croZp(K) = ELX(K)*sideY(K)-ELY(K)*sideX(K)
croMp(K) = croxp(K)*croxp(k)+croyp(K)*croyp(k)+crozp(K)*crozp(k)
cromp(k)=cromp(k)**0.5
c write(6,*)'cros pr',iiii,k,croxp(k),croyp(k),crozp(k),cromp(k)
endif
80117 CONTINUE
if (elx(2).ne.0.and.sidex(1).ne.0) then
croxp(1) = -ELY(2)*sideZ(1)+ELZ(2)*sideY(1)
croyp(1) = +ELX(2)*sideZ(1)-ELZ(2)*sideX(1)
croZp(1) = -ELX(2)*sideY(1)+ELY(2)*sideX(1)
croMp(1) = croxp(1)*croxp(1)+croyp(1)*croyp(1)+crozp(1)*crozp(1)
cromp(1)=cromp(1)**0.5
endif
c write(2,*)'CALCULATION OF SIDEGROUP TORSIONAL ANGLES'
do 8104 k=resno1(iiii)+1,resnum
if (restyp(i).eq.'GLY') goto 8104
if (croxp(k).ne.0.and.crox(k-1).ne.0) then
ara=CROXp(K)*CROX(K-1)+CROYp(K)*CRoy(k-1)+CROZp(K)*CROZ(K-1)
EKFY = -1./CROMp(K)/CROM(K-1)
cosfi = ara*EKFY
if (cosfi.gt.1.) cosfi=1.0
if (cosfi.lt.-1.) cosfi=-1.0
FI = ACOS(cosfi)
gpm = crox(k-1)*sidex(k)+croy(k-1)*sidey(k)+croz(k-1)*sidez(k)
if (gpm.gt.0.0) fi=-fi
faydeg = fi*conv
phip(iiii,k)=faydeg
if (k.eq.10) then
write(6,*)'bond k-1=',elx(k-1),ely(k-1),elz(k-1)
write(6,*)'bond k=',elx(k),ely(k),elz(k)
write(6,*)'bond k+1=',elx(k+1),ely(k+1),elz(k+1)
write(6,*)'side k=',sidex(k),sidey(k),sidez(k)
write(6,*)'crox(k-1)',crox(k-1),croy(k-1),croz(k-1)
write(6,*)'crox(k)',crox(k),croy(k),croz(k)
write(6,*)'croxp(k)',croxp(k),croyp(k),crozp(k)
write(6,*)'bond k+1=',elx(k+1),ely(k+1),elz(k+1)
write(6,*)'torsions:',phi(iiii,k),phip(iiii,k)
write(6,*)'bond angles:',theta(iiii,k),thetpr(iiii,k)
endif
endif
8104 CONTINUE
if (restyp(1).eq.'GLY') goto 8105
if (croxp(1).ne.0.and.crox(2).ne.0) then
ara=CROXp(1)*CROX(2)+CROYp(1)*CRoy(2)+CROZp(1)*CROZ(2)
EKFY = 1./CROMp(1)/CROM(2)
cosfi = ara*EKFY
if (cosfi.gt.1.) cosfi=1.0
if (cosfi.lt.-1.) cosfi=-1.0
FI = ACOS(cosfi)
gpm = crox(2)*sidex(1)+croy(2)*sidey(1)+croz(2)*sidez(1)
if (gpm.gt.0.0) fi=-fi
faydeg = fi*conv
phip(iiii,1)=faydeg
endif
8105 continue
c We also considered terminal side groups. For that purpose we
c assumed that there are two hypothetical backbone bonds, one in
c the same direction as the second virtual bond (elx(3),ely(3),etc)
c such that the first virtual bond has a virtual torsion of 0 deg,
c and the other same direction as bond n-1, such that the terminal
c bond is trans. These are called elx-y-z(1) and elx-y-z(kf),
c kf=resnum+1, and defined in the section following loop 50571.
c***********************************************************************
if(iiii.eq.2) then
do 48933 i=1,resno(iiii)
ax(iiii,i,1)=x(i,2)
ay(iiii,i,1)=y(i,2)
az(iiii,i,1)=z(i,2)
ax(iiii,i,2)=xx(i)
ay(iiii,i,2)=yy(i)
48933 az(iiii,i,2)=zz(i)
endif
open(7,file='cor1cho1.out1')
C COLUMNS: 1 = residue type
c 2-4: x-y-z coordinates of alpha carbons
c 5-7: x-y-z coordinates of side group centroids
write(7,*)resno(iiii)-5
c write(7,*)'backbone and side group coordinates of ',fil(m)
do 20000 i=6,resno(iiii)
20000 write(7,78256)iir(iiii,i),x(i,2),y(i,2),z(i,2),
:xx(i),yy(i),zz(i)
4293 continue
close(7)
c______________________________________________________________________
open(7,file='fi1cho1.out1')
C COLUMNS: 1=residue type
c 2=torsional angle of preceding bond, 3=bond angle,
c 4=side group torsion expressed by preceding bond,
c 5=side group bond length, 6=side group bond angle
do 27990 m=nfile1,nfile2
write(7,*)resno(m)
do 27000 i=7,resno(m)
if(cyss(m,i).eq.1.) iir(m,i)=21
27000 write(7,78255)iir(m,i),phi(m,i),elm(i),theta(m,i)
:,phip(m,i),
:blength(m,i),thetpr(m,i)
27990 continue
close(7)
78255 format(i5,6f9.2)
78256 format(i5,6f9.2)
7824 format(3i5,4f13.4)
c______________________________________________________________________
stop
end