!!this program to add hydrogens into protein structure !!it also can add side chain heavy atoms !13 Jan 2009 intergrate the parameter into this code program main implicit none integer,PARAMETER::NMAX=50000 !!max atom number integer,PARAMETER::Lmax=3000 !!max residue number integer,parameter::NSITT=10 !!maximum clash on the same sites REAL,PARAMETER::pi=3.1415926,hpi=0.017453 !pi and pi/180.0 real,parameter::dismin=1e-6,dismin2=1e-3 !!the minimum cutoff real,DIMENSION(0:NMAX)::xx,yy,zz,xx0,yy0,zz0 real,DIMENSION(NMAX,4)::BOND !record con1E,con1R ... con4E,con4R character*1,dimension(NMAX)::ATY !!record the atom type, H, N O,S, C CHARACTER*4,DIMENSION(NMAX)::ATOMTY,atom !RECORD INPUT ATOMNAME character*3,dimension(NMAX)::resd,resd2 !record the residue sequence name character*4,dimension(20)::afn !full name of 20 amino acids character*30,dimension(NMAX)::fu !for I/O out put integer,dimension(20)::hra,ra !record number of heavy atoms and atoms in each residue !!3, 0 for no form, >0.5 to record on ip(onyl useful on H) of the acceptor integer,DIMENSION(NMAX,5)::bonda !atom connect,5 heavy atom bonded num integer,dimension(NMAX)::ress,ress2 !record residue series number integer,dimension(Lmax,2)::ca !first record Ca atom index, 2nd record residue type in 20 integer,dimension(NMAX)::ASSrd !!!record which atoms are added by this program integer,dimension(Nmax)::tkfp,tkfpres !!for tranfer character*4,dimension(20,24)::atomname integer,dimension(20,24)::atomct !!connect type define in the array character*4,dimension(33)::atomcontype integer,dimension(Nmax)::atc !1record the atomcontype for each atom real,dimension(33,33)::bd0 !!bond length real,dimension(33,33,33)::ag0 !angle integer,dimension(24,4)::bdc integer,dimension(50)::HH1,indx !!this record the atoms index for 1H within 5A, indx the sorted index real,dimension(50)::rHH1 !!this for sorting index CHARACTER*500 SS,du,duu,f1,f2 CHARACTER*54 kfp(Nmax) character*4 abc,prtype character bb real shiftx,shifty,shiftz real dx,dy,dz real rk real vr0(3),vr1(3,3) real ax,ay,az,xi,yi,zi,x,y,z,r real angg,fbd,bd,r1 real r12,r23,r13,bx,by,bz,cx,cy,cz real r2,r3,phi,tab,tac,tbc,t1,ctbcx,ctbcy,ctbcz,ctbax,ctbay,ctbaz,t2 integer ipp(3),ipp0(3) integer ires,iatom,ip1,in,iw,icc,ip2,ip3,ip4 integer ix,iy,iz,ix0,iy0,iz0,lxmin,lymin,lzmin,lxmax,lymax,lzmax,ip integer ia1,ia2,ia3,ia4 integer ia,ib,ic,iflag,ir1,ir2 integer natom0 !!only heavy atoms integer i,j,k,ii,kx,kk,kx0,j1 integer Nfbc,Ledge !the side length of bondy simulation lattice and the edge,smallest box for one residue integer NMAXx,NMAXy,NMAXz !the largest box used in simulation, free boundary condition used:impetrable integer nres,natom,nang,MHbond !number of atoms, angles, hbond bin width integer ires0,ires1 !!!local residue assign and energy calculation use integer iseed,ins,nnzr integer i1,i2,i3,i4,i5,i6,i7,i8 integer iresin,Nhy integer MH1 !the symbol +/- for the H redirection data afn/'GLY','ALA','SER','CYS','VAL','THR','ILE','PRO','MET','ASP','ASN','LEU','LYS','GLU','GLN','ARG','HIS','PHE','TYR','TRP'/ data hra/4,5,6,6,7,7,8,7,8,8,8,8,9,9,9,11,10,11,12,14/ !!heavy atoms data ra/7,10,11,11,16,14,19,14,17,12,14,19,22,15,17,24,17,20,21,24/ !!whole atoms data (atomname(1,i),i=1,7)/" N "," CA "," C "," O "," H ","1HA ","2HA "/ data (atomname(2,i),i=1,10)/" N "," CA "," C "," O "," CB "," H "," HA ","1HB ","2HB ","3HB "/ !ALA data (atomname(3,i),i=1,11)/" N "," CA "," C "," O "," CB "," OG "," H "," HA ","1HB ","2HB "," HG "/ data (atomname(4,i),i=1,11)/" N "," CA "," C "," O "," CB "," SG "," H "," HA ","1HB ","2HB "," HG "/ data (atomname(5,i),i=1,16)/" N "," CA "," C "," O "," CB "," CG1"," CG2"," H "," HA "," HB ","1HG1","2HG1","3HG1","1HG2","2HG2","3HG2"/ data (atomname(6,i),i=1,14)/" N "," CA "," C "," O "," CB "," OG1"," CG2"," H "," HA "," HB "," HG1","1HG2","2HG2","3HG2"/ data (atomname(7,i),i=1,19)/" N "," CA "," C "," O "," CB "," CG1"," CG2"," CD1"," H "," HA "," HB ","1HG1","2HG1","1HG2","2HG2","3HG2","1HD1","2HD1","3HD1"/ data (atomname(8,i),i=1,14)/" N "," CA "," C "," O "," CB "," CG "," CD "," HA ","1HB ","2HB ","1HG ","2HG ","1HD ","2HD "/ data (atomname(9,i),i=1,17)/" N "," CA "," C "," O "," CB "," CG "," SD "," CE "," H "," HA ","1HB ","2HB ","1HG ","2HG ","1HE ","2HE ","3HE "/ !!MET data (atomname(10,i),i=1,12)/" N "," CA "," C "," O "," CB "," CG "," OD1"," OD2"," H "," HA ","1HB ","2HB "/ !!ASP data (atomname(11,i),i=1,14)/" N "," CA "," C "," O "," CB "," CG "," OD1"," ND2"," H "," HA ","1HB ","2HB ","1HD2","2HD2"/ !!ASN data (atomname(12,i),i=1,19)/" N "," CA "," C "," O "," CB "," CG "," CD1"," CD2"," H "," HA ","1HB ","2HB "," HG ","1HD1","2HD1","3HD1","1HD2","2HD2","3HD2"/ !!LEU data (atomname(13,i),i=1,22)/" N "," CA "," C "," O "," CB "," CG "," CD "," CE "," NZ "," H "," HA ","1HB ","2HB ","1HG ","2HG ","1HD ","2HD ","1HE ","2HE ","1HZ ","2HZ ","3HZ "/ !!LYS data (atomname(14,i),i=1,15)/" N "," CA "," C "," O "," CB "," CG "," CD "," OE1"," OE2"," H "," HA ","1HB ","2HB ","1HG ","2HG "/ !!GLU data (atomname(15,i),i=1,17)/" N "," CA "," C "," O "," CB "," CG "," CD "," OE1"," NE2"," H "," HA ","1HB ","2HB ","1HG ","2HG ","1HE2","2HE2"/ !!GLN data (atomname(16,i),i=1,24)/" N "," CA "," C "," O "," CB "," CG "," CD "," NE "," CZ "," NH1"," NH2"," H "," HA ","1HB ","2HB ","1HG ","2HG ","1HD ","2HD "," HE ","1HH1","2HH1","1HH2","2HH2"/ !!ARG data (atomname(17,i),i=1,17)/" N "," CA "," C "," O "," CB "," CG "," ND1"," CD2"," CE1"," NE2"," H "," HA ","1HB ","2HB "," HD1"," HD2"," HE1"/ !!HIS data (atomname(18,i),i=1,20)/" N "," CA "," C "," O "," CB "," CG "," CD1"," CD2"," CE1"," CE2"," CZ "," H "," HA ","1HB ","2HB "," HD1"," HD2"," HE1"," HE2"," HZ "/ !!PHE data (atomname(19,i),i=1,21)/" N "," CA "," C "," O "," CB "," CG "," CD1"," CD2"," CE1"," CE2"," CZ "," OH "," H "," HA ","1HB ","2HB "," HD1"," HD2"," HE1"," HE2"," HH "/ !!TYR data (atomname(20,i),i=1,24)/" N "," CA "," C "," O "," CB "," CG "," CD1"," CD2"," NE1"," CE2"," CE3"," CZ2"," CZ3"," CH2"," H "," HA ","1HB ","2HB "," HD1"," HE1"," HE3"," HZ2"," HZ3"," HH2"/ !!TRP !!assign the atom connection type total 29 types data (atomct(1,i),i=1,7)/10,2,4,26,17,19,19/ data (atomct(2,i),i=1,10)/10,1,4,26,3,17,19,18,18,18/ !ALA data (atomct(3,i),i=1,11)/10,1,4,26,2,28,17,19,18,18,17/ !SER data (atomct(4,i),i=1,11)/10,1,4,26,2,29,17,19,18,18,21/ !CYS data (atomct(5,i),i=1,16)/10,1,4,26,1,3,3,17,19,18,18,18,18,18,18,18/ !VAL data (atomct(6,i),i=1,14)/10,1,4,26,1,28,3,17,19,18,17,18,18,18/ !THR data (atomct(7,i),i=1,19)/10,1,4,26,1,2,3,3,17,19,18,18,18,18,18,18,18,18,18/ !ILE data (atomct(8,i),i=1,14)/9,30,4,26,31,31,32,19,18,18,18,18,18,18/ !PRO data (atomct(9,i),i=1,17)/10,1,4,26,2,2,29,3,17,19,18,18,18,18,18,18,18/ !!MET data (atomct(10,i),i=1,12)/10,1,4,26,2,5,27,27,17,19,18,18/ !!ASP data (atomct(11,i),i=1,14)/10,1,4,26,2,5,26,11,17,19,18,18,17,17/ !!ASN data (atomct(12,i),i=1,19)/10,1,4,26,2,1,3,3,17,19,18,18,18,18,18,18,18,18,18/ !!LEU data (atomct(13,i),i=1,22)/10,1,4,26,2,2,2,2,12,17,19,18,18,18,18,18,18,18,18,20,20,20/ !!LYS data (atomct(14,i),i=1,15)/10,1,4,26,2,2,5,27,27,17,19,18,18,18,18/ !!GLU data (atomct(15,i),i=1,17)/10,1,4,26,2,2,5,26,11,17,19,18,18,18,18,17,17/ !!GLN data (atomct(16,i),i=1,24)/10,1,4,26,2,2,2,13,4,13,13,17,19,18,18,18,18,18,18,20,20,20,20,20/ !!ARG data (atomct(17,i),i=1,17)/10,1,4,26,2,6,15,6,7,15,17,19,18,18,17,22,22/ !!HIS data (atomct(18,i),i=1,20)/10,1,4,26,2,23,23,23,23,23,23,17,19,18,18,25,25,25,25,25/ !!PHE data (atomct(19,i),i=1,21)/10,1,4,26,2,23,23,23,23,23,23,28,17,19,18,18,25,25,25,25,17/ !!TYR data (atomct(20,i),i=1,24)/10,1,4,26,2,14,23,8,16,8,23,23,23,23,17,19,18,18,25,17,25,25,25,25/ !!TRP !! 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 data atomcontype/"CT1","CT2","CT3","C","CC","CPH1","CPH2","CPT","N","NH1","NH2","NH3","NC2","CY","NR2","NY","H","HA","HB","HC","HS","HR1","CA","HR3","HP","O","OC","OH1","S","CP1","CP2","CP3","NP"/ ! atom types 8+8+9+4 !18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33(P_Nterm) !!bond length assign in 29*29 matrix, angle anssign in 29*29*29 matrix f1='' call getarg(1,f1) if(len(trim(f1))<2)then print*,"USING Guide::" print*,"./haad xx.pdb (ADD H only, output: xx.pdb.h)" goto 777 endif !!read pdb structure *************************FULL ATOMIC*************************************** 1101 format(A100) 1102 format(A30) open(1,file=trim(f1)) !!!only for test in=0 ir1=1 30 read(1,1101,end=31)ss !!F1 structure file if(ss(1:3).eq.'END'.or.ss(1:3).eq.'TER') goto 31 bb=ss(17:17) if(ss(1:4).eq.'ATOM'.and.(bb.eq.' '.or.bb.eq.'A'))then in=in+1 read(ss,101)du,atom(in),duu,resd(in),duu,ir2,duu,xx0(in),yy0(in),zz0(in) read(ss,1102)kfp(in) tkfpres(in)=ir2 if(ir2>ic.and.in>1)ir1=ir1+1 ress(in)=ir1 resd2(ir1)=resd(in) !!get the residue sequence ic=ir2 endif goto 30 31 continue close(1) if(in<10)then print*,"read pdb file err, total atoms",in stop endif natom0=in Nres=ir1 !!!according to ress build the force field !!assign the known coordinates into xyz call getresin !assign the atomname, atom index and CA xx=999.0 yy=999.0 zz=999.0 tkfp=0 do i=1,natom !!should be charmm atm file do j=1,natom0 if(atomty(i).eq.atom(j).and.ress2(i).eq.ress(j))then xx(i)=xx0(j) yy(i)=yy0(j) zz(i)=zz0(j) tkfp(i)=j goto 11 endif enddo 11 continue enddo !!count the number of heavy atoms bonded to each atom do i=1,natom kx0=0 do k=1,4 ii=bonda(i,k) if(ii>0)then if(aty(ii).eq.'H')cycle kx0=kx0+1 endif enddo bonda(i,5)=kx0 !number of heavy atoms bonded ress(i)=ress2(i) enddo !do i=1,Natom ! write(*,"A4,5I5,A4,A5"),"tstk",i,ress(i),ress2(i),ca(ress(i),1),ca(ress(i),2),resd(i),atomty(i) !enddo Assrd=0 iseed=natom 101 format(A12,A4,A1,A3,A2,I4,A4,3F8.3) !pdb 772 continue call heavy_ADD call HADD f2=trim(f1)//".h" open(10,file=f2) do i=1,natom k=tkfp(i) ! if(xx(i)>990.0)cycle if(k>1)ir2=tkfpres(k) if(ASSrd(i)>0)then write(10,129),0,atomty(i),resd(i),ir2,xx(i),yy(i),zz(i),ASSrd(i) elseif(k>0)then write(10,1103),kfp(k),xx(i),yy(i),zz(i) endif enddo write(10,119) close(10) 1103 format(A30,3f8.3) 119 format("TER") 129 format("ATOM",I7,1X,A4,1X,A3,2X,I4,4X,3F8.3,I3) !pdb 777 continue contains !!6 and 66 done !!this subroutine to add those misssing heavy atoms !!except CB, CG and CD subroutine HEAVY_ADD implicit none !scan heavy atoms at first do i=1,natom if(aty(i).eq.'H'.or.xx(i)<990)cycle ASSrd(i)=1 ires=ress(i) do k=1,4 ii=bonda(i,k) if(ii<1)goto 51 if(xx(ii)>990)cycle icc=ii bd=bond(i,k) goto 51 enddo 51 continue ax=xx(icc) !!origin point ay=yy(icc) az=zz(icc) !write(*,"A10,I5,A4,A4,I4,3f8.3"),"readd heavy",i,atomty(i),resd(i),ress(i),ax,ay,az kx0=0 ip2=0 do j=1,4 !!the connect heavy atom !!ipx-ip2-icc-H (1,2,3) ii=bonda(icc,j) if(ii<1)goto 301 if(xx(ii)>990.or.ii.eq.i.or.aty(ii).eq."H")cycle ip2=ii kx0=kx0+1 ipp0(kx0)=ip2 enddo 301 continue i1=atc(i) i2=atc(icc) i3=atc(ip2) if(i1>i3)then i1=atc(ip2) i3=atc(i) endif angg=ag0(i1,i2,i3)*hpi kx=0 do j=1,4 ii=bonda(ip2,j) if(ii<1.or.(ii==icc.or.aty(ii)=='H'.or.xx(ii)>990))cycle kx=kx+1 ipp(kx)=ii enddo ip3=ipp(1) !treat heavy atoms missing !ip3-ip2-icc-i (heavy atom) !!C-term ALA-Cb; ILE Cd; LEU Cd1,Cd2;MET Ce; Val CG1; THR CG2 sp3 !!O-term O for all, ASN OD1; ASP OD1, OD2; GLN OE1; GLU OE1,OE2; TYR OH;sp2 || SER OG; THR OG1, OXT (sp3) !!S-term CYS SG sp3 !!N-term N-term LYS NZ sp3, ARG NH1,NH2; ASN ND2; GLN NE2 sp2 abc=atomty(i) !for sp2 second atom if(bonda(i,5).eq.1)then if(abc.eq.' OD2'.or.abc.eq.' OE2'.or.abc.eq.' ND2'.or.abc.eq.' NE2'.or.abc.eq.' NH2'.or.abc.eq.' OH '.or.abc.eq.' OXT')then !print*,"addheavy",abc,i if(abc.eq.' OH ')then ip2=icc-1 ip3=icc-2 else !if(abc.eq." OXT".or.abc.eq." NH2")then ip2=icc-1 ip3=icc+1 endif x=0.5*(xx(ip2)+xx(ip3)) y=0.5*(yy(ip2)+yy(ip3)) z=0.5*(zz(ip2)+zz(ip3)) r1=sqrt((x-ax)*(x-ax)+(y-ay)*(y-ay)+(z-az)*(z-az)) if(r1>dismin2)then fbd=bd/r1 xx(i)=ax+fbd*(ax-x) yy(i)=ay+fbd*(ay-y) zz(i)=az+fbd*(az-z) else ASSrd(i)=2 ! print*,"warn!ADD Cterm1",i,atomty(i),ip3,ip2,icc ip3=ipp(1) xx(i)=ax+xx(ip2)-xx(ip3) yy(i)=ay+yy(ip2)-yy(ip3) zz(i)=az+zz(ip2)-zz(ip3) endif elseif(abc.eq.' CD2'.and.resd(i).eq.'LEU')then ! LEU xi=ax-xx(icc-1) yi=ay-yy(icc-1) zi=az-zz(icc-1) x=xx(i-1) !CD1 y=yy(i-1) z=zz(i-1) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,2.18166) !125o xx(i)=x yy(i)=y zz(i)=z elseif(abc.eq.' CG2')then !VAL, THR, ILE xi=ax-xx(icc-3) yi=ay-yy(icc-3) zi=az-zz(icc-3) x=xx(i-1) !CD1 y=yy(i-1) z=zz(i-1) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,2.18166) !125o xx(i)=x yy(i)=y zz(i)=z ! print*,"CG2",atomty(i),xx(i),yy(i),zz(i) elseif((abc.eq.' CD1'.and.resd(i).eq."ILE").or.(abc.eq.' CE '.and.resd(i).eq."MET"))then !ILE xi=xx(icc-1)-0.5*(xx(i-1)+xx(icc-4)) yi=yy(icc-1)-0.5*(yy(i-1)+yy(icc-4)) zi=zz(icc-1)-0.5*(zz(i-1)+zz(icc-4)) x=xi-ax y=yi-ay z=zi-az r1=sqrt(x*x+y*y+z*z) if(r1<0.1)r1=0.1 fbd=bd/r1 xx(i)=ax+fbd*x yy(i)=ay+fbd*y zz(i)=az+fbd*z else !!sp3 and sp2 first atom ????????? if(abc.eq.' CG1'.or.abc.eq.' SG '.or.abc.eq.' OG1'.or.abc.eq.' OG ')ip3=ip2+1 ! print*,"HADD>=",ress(i),atomty(ip3),atomty(ip2),atomty(i),resd(i),abc x=xx(ip2)-xx(ip3) y=yy(ip2)-yy(ip3) z=zz(ip2)-zz(ip3) r1=sqrt(x*x+y*y+z*z) if(r1>dismin2)then fbd=bd/r1 xx(i)=ax+fbd*x yy(i)=ay+fbd*y zz(i)=az+fbd*z else ASSrd(i)=2 ! print*,"warn!ADD Cterm2",i,atomty(i),ip3,ip2,icc xx(i)=2.0*ax-xx(ip2) yy(i)=2.0*ay-yy(ip2) zz(i)=2.0*az-zz(ip2) endif endif elseif(bonda(i,5)>=2)then !!treat on more atoms missing if(abc.eq.' CG ')ip3=ip2+1 x=xx(ip2)-xx(ip3) y=yy(ip2)-yy(ip3) z=zz(ip2)-zz(ip3) r1=sqrt(x*x+y*y+z*z) if(r1>dismin2)then fbd=bd/r1 xx(i)=ax+fbd*x yy(i)=ay+fbd*y zz(i)=az+fbd*z else ASSrd(i)=2 ! print*,"warn!ADD Cterm3",i,atomty(i),ip3,ip2,icc xx(i)=ax+x yy(i)=ay+y zz(i)=az+z endif !!default set if(.not.(ca(ress(i),2).eq.8.or.ca(ress(i),2)>=17))cycle !print*,atomty(i),resd(i),ress(i),"<<<dismin2)then fbd=bd/r1 xx(i)=ax+fbd*(ax-x) yy(i)=ay+fbd*(ay-y) zz(i)=az+fbd*(az-z) else ASSrd(i)=2 ip3=ipp(1) xx(i)=ax+xx(ip2)-xx(ip3) yy(i)=ay+yy(ip2)-yy(ip3) zz(i)=az+zz(ip2)-zz(ip3) endif endif endif !!ring finish endif end if !print*,i,atomty(i),resd(i),xx(i),yy(i),zz(i) !if(xx(i)>990.0)print*,i,atomty(i),resd(i),xx(i) enddo !!heavy atom scan finish !print*,"HEAVY add finish" end subroutine HEAVY_ADD !!this subroutine add hydrogen atoms subroutine HADD implicit none !call HEAVY_ADD i=natom+1 do while(i>2) i=i-1 if(xx(i)<990.or.aty(i)/='H')cycle icc=bonda(i,1) !!bonded heavy atom ax=xx(icc) !!origin point ay=yy(icc) az=zz(icc) bd=bond(i,1) ires=ress(i) kx0=0 ! bd=bd*0.864 do j=1,4 !!the connect heavy atom !!ipx-ip2-icc-H (1,2,3) ii=bonda(icc,j) if(ii<1)goto 301 if(xx(ii)>990.or.ii.eq.i.or.aty(ii).eq."H")cycle ip2=ii kx0=kx0+1 ipp0(kx0)=ip2 enddo 301 continue ip2=ipp0(1) ASSrd(i)=1 !!!get the angle ip2-icc-i i1=atc(i) i2=atc(icc) i3=atc(ip2) if(i1>i3)then i1=atc(ip2) i3=atc(i) endif angg=ag0(i1,i2,i3)*hpi kx=0 do j=1,4 ii=bonda(ip2,j) if(ii<1.or.ii==icc.or.ATY(ii)=='H'.or.xx(ii)>990)cycle ip3=bonda(ip2,j) kx=kx+1 ipp(kx)=ip3 enddo vr0(1)=xx(ip2)-ax vr0(2)=yy(ip2)-ay vr0(3)=zz(ip2)-az !if(ca(ires,2).eq.1)then ! write(*,"I5,2I4,A4,A5,3I5,2f6.3"),i,ires,ca(ires,2),resd(i),atomty(i),icc,ip2,ip3,bd,angg !endif !!!rotation axis , cross of vr0 and vr1, first H is transto ip2 if(atomty(i)(1:2)=='3H')then !!term sp3-sp3 trans-, Gauche-, Gauche+ if(aty(icc).eq."C")bd=bd*0.933 !1.05 ASSrd(i-1)=1 ASSrd(i-2)=1 do k=1,kx vr1(k,1)=xx(ipp(k))-xx(ip2) !!connect heavy atoms vectors vr1(k,2)=yy(ipp(k))-yy(ip2) !!connect heavy atoms vectors vr1(k,3)=zz(ipp(k))-zz(ip2) !!connect heavy atoms vectors enddo x=xx(ip2) y=yy(ip2) z=zz(ip2) xi=vr1(1,2)*vr0(3)-vr1(1,3)*vr0(2) yi=vr1(1,3)*vr0(1)-vr1(1,1)*vr0(3) zi=vr1(1,1)*vr0(2)-vr1(1,2)*vr0(1) !print*,"ccc",x,y,z,xi,yi,zi,ax,ay,az,angg call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,angg) !H on N, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>5.0)then ! print*, "warning CH3-NH3-1",r1,bd,ress2(i),i,atom(i,1) r1=sqrt((xx(ip3)-xx(ip2))**2+(yy(ip3)-yy(ip2))**2+(zz(ip3)-zz(ip2))**2) if(r1>10.0.or.r1<0.001)then ! print*,"neighbour position severa wrong CH3-NH3, again" xx(i)=ax+(xx(ip2)-xx(ip3))/bond(ip2,2) yy(i)=ay+(yy(ip2)-yy(ip3))/bond(ip2,2) zz(i)=az+(zz(ip2)-zz(ip3))/bond(ip2,2) ASSrd(i)=2 else fbd=bd/r1 xx(i)=ax+(xx(ip2)-xx(ip3))*fbd yy(i)=ay+(yy(ip2)-yy(ip3))*fbd zz(i)=az+(zz(ip2)-zz(ip3))*fbd endif else !print*,"ddd",x,y,z,xi,yi,zi,ax,ay,az,r1,bd,angg fbd=bd/r1 !!adjust bond length xx(i)=ax+fbd*(x-ax) yy(i)=ay+fbd*(y-ay) zz(i)=az+fbd*(z-az) endif angg=2.0*pi/3.0 x=xx(i) y=yy(i) z=zz(i) xi=vr0(1) yi=vr0(2) zi=vr0(3) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,angg) !H on N, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>100.0)then ! print*, "warning CH3,NH3-2",r1,bd,ress2(i),i,atom(i,1) xx(i-1)=ax+cos(109.5/180*pi)*(xx(i)-ax) yy(i-1)=ay+sin(109.5/180*pi)*(yy(i)-ay) zz(i-1)=az ASSrd(i-1)=2 else fbd=bd/r1 !!adjust bond length xx(i-1)=ax+fbd*(x-ax) yy(i-1)=ay+fbd*(y-ay) zz(i-1)=az+fbd*(z-az) endif x=xx(i) y=yy(i) z=zz(i) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,-angg) !H on N, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>10.0)then ! print*, "warning CH3,NH3-2",r1,bd,ress2(i),i,atom(i,1) xx(i-2)=ax+cos(109.5/180*pi)*(xx(i)-ax) yy(i-2)=ay+sin(-109.5/180*pi)*(yy(i)-ay) zz(i-2)=az ASSrd(i-1)=2 else fbd=bd/r1 !!adjust bond length xx(i-2)=ax+fbd*(x-ax) yy(i-2)=ay+fbd*(y-ay) zz(i-2)=az+fbd*(z-az) endif i=i-2 !!!!sp3 -term finish **************************************************************** elseif(atomty(i)(1:2)=='2H')then !!in term or middle ! bd=bd*1.08 if(aty(icc).eq."C")bd=bd*0.933 ASSrd(i-1)=1 if(ATY(icc)=='N')then !!in term, sp2 ! if(kx==2)then !!sp2-sp2 term if(ATY(ipp(1)).eq.'C')then !!N=O -C ip3=ipp(2) else ip3=ipp(1) endif !!the O ! if(resd(i).eq."ARG")then !!N-C-N-H ! angg=1.02*angg ! else !!O=C-N-H ! angg=0.98*angg ! endif vr1(1,1)=xx(ip3)-xx(ip2) vr1(1,2)=yy(ip3)-yy(ip2) vr1(1,3)=zz(ip3)-zz(ip2) x=xx(ip2) y=yy(ip2) z=zz(ip2) xi=vr1(1,2)*vr0(3)-vr1(1,3)*vr0(2) yi=vr1(1,3)*vr0(1)-vr1(1,1)*vr0(3) zi=vr1(1,1)*vr0(2)-vr1(1,2)*vr0(1) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,angg) !H on N, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>5.0)then ASSrd(i)=2 r1=sqrt((xx(ip3)-xx(ip2))**2+(yy(ip3)-yy(ip2))**2+(zz(ip3)-zz(ip2))**2) if(r1>10.0.or.r1<0.001)then ! print*,"neighbour bond position severa wrong NH2 again" xx(i)=ax+(xx(ip2)-xx(ip3))*0.6 yy(i)=ay+(yy(ip2)-yy(ip3))*0.6 zz(i)=az+(zz(ip2)-zz(ip3))*0.6 else fbd=bd/r1 xx(i)=ax+(xx(ip2)-xx(ip3))*fbd yy(i)=ay+(yy(ip2)-yy(ip3))*fbd zz(i)=az+(zz(ip2)-zz(ip3))*fbd endif else fbd=bd/r1 !!adjust bond length xx(i)=ax+fbd*(x-ax) yy(i)=ay+fbd*(y-ay) zz(i)=az+fbd*(z-az) endif x=xx(ip2) y=yy(ip2) z=zz(ip2) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,-angg) !H on N, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>5.0)then r1=sqrt((xx(ip3)-xx(ip2))**2+(yy(ip3)-yy(ip2))**2+(zz(ip3)-zz(ip2))**2) if(r1>5.0.or.r1<0.001)then xx(i-1)=ax+(xx(ip2)-xx(ip3))*0.6 yy(i-1)=ay+(yy(ip2)-yy(ip3))*0.6 zz(i-1)=az+(zz(ip2)-zz(ip3))*0.6 ASSrd(i-1)=2 print*,"severa wrong NH2-2",i-1,xx(i-1),yy(i-1) else fbd=bd/r1 xx(i-1)=ax+(xx(ip3)-xx(ip2))*fbd yy(i-1)=ay+(yy(ip3)-yy(ip2))*fbd zz(i-1)=az+(zz(ip3)-zz(ip2))*fbd endif else fbd=bd/r1 !!adjust bond length xx(i-1)=ax+fbd*(x-ax) yy(i-1)=ay+fbd*(y-ay) zz(i-1)=az+fbd*(z-az) endif if(ca(i,2)==8.and.ress2(i)==1)then !!PRO N-term xi=vr0(1) yi=vr0(2) zi=vr0(3) do k=-1,0 x=xx(i+k) y=yy(i+k) z=zz(i+k) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,0.5*pi) !H on N, each H is trans xx(i+k)=x yy(i+k)=y zz(i+k)=z enddo endif else !!!C inter connected, not in term, rotate first H on of ip2 around icc-ip4 if(ipp0(2)==ip2)then ip4=ipp0(1) if(abs(icc-ip4)>20)print*,"ip4 wrong!!!" else ip4=ipp0(2) endif angg=1.1*angg !!!120/110 degree do ii=1,2 if(ii==2)then i2=ip2 ip2=ip4 ip4=i2 angg=-angg endif x=xx(ip2) y=yy(ip2) z=zz(ip2) xi=xx(ip4)-ax yi=yy(ip4)-ay zi=zz(ip4)-az call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,angg) !H on C, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>5.0)then r1=sqrt((xx(ip3)-xx(ip2))**2+(yy(ip3)-yy(ip2))**2+(zz(ip3)-zz(ip2))**2) if(r1>10.0.or.r1<0.001)then ! print*,"neighbour bond position severa wrong CH2-1 again" xx(i)=ax+(xx(ip2)-xx(ip3))/bond(ip2,2) yy(i)=ay+(yy(ip2)-yy(ip3))/bond(ip2,2) zz(i)=az+(zz(ip2)-zz(ip3))/bond(ip2,2) ASSrd(i)=2 else fbd=bd/r1 xx(i)=ax+(xx(ip3)-xx(ip2))*fbd yy(i)=ay+(yy(ip3)-yy(ip2))*fbd zz(i)=az+(zz(ip3)-zz(ip2))*fbd endif else fbd=bd/r1 !!adjust bond length if(ii<2)then xx(i)=ax+fbd*(x-ax) yy(i)=ay+fbd*(y-ay) zz(i)=az+fbd*(z-az) else xx(i)=0.5*(xx(i)+ax+fbd*(x-ax)) yy(i)=0.5*(yy(i)+ay+fbd*(y-ay)) zz(i)=0.5*(zz(i)+az+fbd*(z-az)) endif endif x=xx(ip2) y=yy(ip2) z=zz(ip2) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,-angg) !H on C, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>5.0)then r1=sqrt((xx(ip3)-xx(ip2))**2+(yy(ip3)-yy(ip2))**2+(zz(ip3)-zz(ip2))**2) if(r1>10.0.or.r1<0.001)then print*,"neighbour bond position severa wrong CH2-2",i,xx(i) xx(i)=ax+(xx(ip2)-xx(ip3))/bond(ip2,2) yy(i)=ay+(yy(ip2)-yy(ip3))/bond(ip2,2) zz(i)=az+(zz(ip2)-zz(ip3))/bond(ip2,2) ASSrd(i)=2 else fbd=bd/r1 xx(i-1)=ax+(xx(ip2)-xx(ip3))*fbd yy(i-1)=ay+(yy(ip2)-yy(ip3))*fbd zz(i-1)=az+(zz(ip2)-zz(ip3))*fbd endif else fbd=bd/r1 !!adjust bond length if(ii<2)then xx(i-1)=ax+fbd*(x-ax) yy(i-1)=ay+fbd*(y-ay) zz(i-1)=az+fbd*(z-az) else xx(i-1)=0.5*(xx(i-1)+ax+fbd*(x-ax)) yy(i-1)=0.5*(yy(i-1)+ay+fbd*(y-ay)) zz(i-1)=0.5*(zz(i-1)+az+fbd*(z-az)) endif endif enddo !!!ii=1,2 endif !! 2H finished i=i-1 else !!1H case, sp2 has been assign, O-H and S-H also done trans-; only for sp3 C case if(kx0==3)then !!sp3 case bd=bd*0.945 x=(xx(ipp0(1))+xx(ipp0(2))+xx(ipp0(3)))/3.0 y=(yy(ipp0(1))+yy(ipp0(2))+yy(ipp0(3)))/3.0 z=(zz(ipp0(1))+zz(ipp0(2))+zz(ipp0(3)))/3.0 r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>10.0)then ! print*, "warning H1-sp3",ress2(i),atom(i,1),i,r1 ASSrd(i)=2 xx(i)=ax+0.3 yy(i)=ay+0.6 zz(i)=az+0.7 else fbd=bd/r1 !!adjust bond length xx(i)=ax+fbd*(ax-x) yy(i)=ay+fbd*(ay-y) zz(i)=az+fbd*(az-z) endif elseif(kx0==2)then !!sp2, aromatic ring H, H-N bd=bd*0.95 x=(xx(ipp0(1))+xx(ipp0(2)))*0.5 y=(yy(ipp0(1))+yy(ipp0(2)))*0.5 z=(zz(ipp0(1))+zz(ipp0(2)))*0.5 r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>5.0)then ASSrd(i)=2 xx(i)=ax+0.3 yy(i)=ay+0.6 zz(i)=az+0.7 else fbd=bd/r1 !!adjust bond length xx(i)=ax+fbd*(ax-x) yy(i)=ay+fbd*(ay-y) zz(i)=az+fbd*(az-z) endif elseif(kx0==1)then !!H-O,H-S bd=bd*0.94 vr1(1,1)=xx(ipp(1))-xx(ip2) vr1(1,2)=yy(ipp(1))-yy(ip2) vr1(1,3)=zz(ipp(1))-zz(ip2) x=xx(ip2) y=yy(ip2) z=zz(ip2) xi=vr1(1,2)*vr0(3)-vr1(1,3)*vr0(2) yi=vr1(1,3)*vr0(1)-vr1(1,1)*vr0(3) zi=vr1(1,1)*vr0(2)-vr1(1,2)*vr0(1) angg=angg call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,angg) !H on N, each H is trans r1=sqrt((x-ax)**2+(y-ay)**2+(z-az)**2) if(r1<0.01.or.r1>5.0)then ! print*, "warning O-H",ress2(i),atom(i,1),i,r1 ASSrd(i)=2 xx(i)=ax+0.3 yy(i)=ay+0.6 zz(i)=az+0.7 else fbd=bd/r1 !!adjust bond length xx(i)=ax+fbd*(x-ax) yy(i)=ay+fbd*(y-ay) zz(i)=az+fbd*(z-az) endif endif endif !!3H enddo !!kk144.0)then if(ir1>",i,j,r,np0 cycle else goto 231 endif endif if(r>0.01.and.r<25.0)then np0=np0+1 HH1(np0)=j rHH1(np0)=r !!use aty(j) to judge whether its overlap endif enddo 231 continue if(np0<1)cycle ! print*,"nop",i,np0,atomty(i) do j=1,np0 arrin(j)=rHH1(j) enddo if(np0>1)then call sortindex(np0,arrin,indx) !!increaing order else indx(1)=1 endif !!decide the vecoter summ x=0 y=0 z=0 do j=1,np0 !!four cases, Hbond yes/not, attractive/repulsive ind=indx(j) ip=HH1(ind) !!atom id if(rHH1(ind)<0.01)cycle if(aty(ip).eq."C".or.aty(ip).eq."S")then r=9.0 !3.0^2 elseif(aty(ip).eq."N")then r=8.12 !2.85^2 elseif(aty(ip).eq."H")then r=4.0 elseif(aty(ip).eq."O")then r=7.67 !2.77^2 endif if(rHH1(ind)=1.0)cycle angg=acos(rf) x=xx(i) y=yy(i) z=zz(i) call rotate_matrix2(x,y,z,ax,ay,az,xi,yi,zi,angg) !H on N, each H is trans ! r=sqrt((x-xx(i))**2+(y-yy(i))**2+(z-zz(i))**2) ! r=sqrt((x-xx(icc))**2+(y-yy(icc))**2+(z-zz(icc))**2) !write(*,"A10,A5,I5,f6.3,7f8.3"),"rotation angg",atomty(i),i,angg,xx(i),yy(i),zz(i),x,y,z,r xx(i)=x yy(i)=y zz(i)=z endif !!H end do !!all spH finish end subroutine AD_1H !this sub for sort, output in increase series ! puts into a the permutation vector which sorts v into increasing order. SUBROUTINE sort(RA,N) implicit none real RA(N) real rra integer N,ir,L,i,j L=N/2+1 IR=N 10 IF (L.GT.1) THEN L=L-1 RRA=RA(L) ELSE RRA=RA(IR) RA(IR)=RA(1) IR=IR-1 IF (IR.EQ.1) THEN RA(1)=RRA RETURN END IF END IF I=L J=L+L 20 IF (J.LE.IR) THEN IF (J.LT.IR) THEN IF (RA(J).LT.RA(J+1)) J=J+1 END IF IF (RRA.LT.RA(J)) THEN RA(I)=RA(J) I=J J=J+J ELSE J=IR+1 END IF GO TO 20 END IF RA(I)=RRA GO TO 10 END subroutine sort !!this subroutine to sort array a and return accerating index after sorting !!Note the input array will change value SUBROUTINE sortindex(n,arrin,indx) implicit none integer n integer indx(n) real arrin(n) integer i,j,k,ir,indxt real q do j=1,n indx(j)=j end do k=n/2+1 ir=n do if(k>1)then k=k-1 indxt=indx(k) q=arrin(indxt) else indxt=indx(ir) q=arrin(indxt) indx(ir)=indx(1) ir=ir-1 if(ir==1) then indx(1)=indxt return endif endif i=k j=k+k do while(j<=ir) if(jNhy)then if(bdc(iatom,j)>3)then bonda(iatom,j)=iatom+bdc(iatom,j)+2 else bonda(iatom,j)=iatom+bdc(iatom,j) endif enddo enddo bonda(1,3)=Nhy+1 bonda(1,4)=Nhy+2 bonda(iatom,1)=1 bonda(iatom+1,1)=1 iatom=iatom+1 do i=Nhy+1,nat iatom=iatom+1 bonda(iatom,1)=iatom+bdc(i,1)-2 enddo bonda(3,2)=ca(2,1)-1 if(Iresin.eq.8)then !!PRO bonda(1,2)=Nhy endif !!!N_term finish do k=2,Nres-1 iresin=ca(k,2) bdc=0 nat=ra(iresin) Nhy=hra(iresin) call Assbonda do i=1,nat iatom=iatom+1 do j=1,4 if(bdc(i,j)==0)cycle bonda(iatom,j)=iatom+bdc(i,j) enddo enddo ipca=ca(k,1) bonda(ipca-1,3)=ca(k-1,1)+1 bonda(ipca+1,2)=ca(k+1,1)-1 enddo !C-Term iresin=ca(Nres,2) bdc=0 nat=ra(iresin) Nhy=hra(iresin) ipca=ca(Nres,1) call Assbonda do i=1,Nhy iatom=iatom+1 do j=1,4 if(bdc(i,j)==0)cycle if(bdc(i,j)>3)then bonda(iatom,j)=iatom+bdc(i,j)+1 else bonda(iatom,j)=iatom+bdc(i,j) endif enddo enddo bonda(iatom+1,1)=ipca+1 bonda(ipca+1,2)=iatom+1 iatom=iatom+1 do i=Nhy+1,nat iatom=iatom+1 bonda(iatom,1)=iatom+bdc(i,1)-1 enddo bonda(ipca-1,3)=ca(Nres-1,1)+1 if(iresin.eq.8)then bonda(ca(Nres,1)-1,2)=bonda(ca(Nres,1)-1,2)-1 endif 333 continue !bond vector assign finish do i=1,0 !Natom nat=bonda(i,2) if(nat>=1)then write(*,1104),i,atomty(i),resd(i),bonda(i,1),bonda(i,2),bonda(i,3),bonda(i,4),atomty(bonda(i,1)),atomty(nat),atomty(bonda(i,3)),atomty(bonda(i,4)) else write(*,1105),i,atomty(i),resd(i),bonda(i,1),bonda(i,2),bonda(i,3),bonda(i,4),atomty(bonda(i,1)),nat endif enddo 1104 format(I4,A5,A4,4I5,4A5) 1105 format(I4,A5,A4,4I5,A5,I3) !!assign the bondlength call bondp do i=1,natom i1=atc(i) do j=1,4 if(bonda(i,j)<1)cycle i2=atc(bonda(i,j)) if(i1<=i2)then bond(i,j)=bd0(i1,i2) else bond(i,j)=bd0(i2,i1) endif if(bond(i,j)<0.5)then write(*,1106),"bondm",i,bonda(i,j),i1,i2,atomty(i),atomty(bonda(i,j)),bd0(i1,i2),bd0(i2,i1) endif enddo enddo 1106 format(A6,4I5,2A5,2f8.3) !!assign angle ///********************////////////// call angp do i=1,natom !!check whether all are in position ! if(aty(i).ne."H")cycle kx0=0 ip2=0 icc=bonda(i,1) do j=1,4 !!the connect heavy atom !!ipx-ip2-icc-H (1,2,3) ii=bonda(icc,j) if(ii<1)goto 301 if(xx(ii)>990.or.ii.eq.i.or.aty(ii).eq."H")cycle ip2=ii kx0=kx0+1 ipp0(kx0)=ip2 enddo 301 continue if(ip2<1)goto 302 i1=atc(i) i2=atc(icc) i3=atc(ip2) if(i1>i3)then i1=atc(ip2) i3=atc(i) endif angg=ag0(i1,i2,i3) if(angg<1.0)then write(*,1107),"anm",i,icc,ip2,atomty(i),atomty(icc),atomty(ip2),i1,i2,i3,atomcontype(i1),atomcontype(i2),atomcontype(i3),bonda(icc,1),bonda(icc,2),bonda(icc,3),bonda(icc,4) endif 302 continue enddo 1107 format(A6,3I5,3A5,3I3,3A5,4I5) end subroutine getresin !!this subroutine assign the bond network !!which atoms are bonded subroutine assbonda implicit none integer ig,nratom !!residue type !!seperate backbone and sidechain !!the internal difference Nhy=hra(iresin) !print*,iresin,Nhy bdc(1,1)=1 bdc(1,2)=Nhy !1,3 to C bdc(2,1)=3 bdc(2,2)=-1 bdc(2,3)=1 bdc(2,4)=Nhy bdc(3,1)=-1 bdc(3,3)=1 ! bdc(3,4)=1 !bdc4 to N !!for backbone, only GLY and PRO special bdc(4,1)=-1 ! bdc(4,2)=-1 bdc(Nhy+1,1)=-Nhy !H bdc(Nhy+2,1)=-Nhy !HA bdc(5,1)=-3 !CB-Ca bdc(5,2)=1 bdc(5,3)=Nhy-1 bdc(5,4)=Nhy-2 bdc(Nhy+3,1)=-Nhy+2 bdc(Nhy+4,1)=-Nhy+1 if(iresin.eq.1)then !GLY bdc(2,1)=Nhy bdc(2,4)=Nhy+1 bdc(Nhy+1,1)=-Nhy bdc(5,2)=0 bdc(5,3)=0 bdc(5,4)=0 bdc(nhy+2,1)=-Nhy bdc(nhy+3,1)=-Nhy-1 elseif(iresin.eq.2)then !ALA bdc(Nhy,2)=3 bdc(Nhy,3)=4 bdc(Nhy,4)=5 bdc(Nhy+3,1)=-Nhy+2 bdc(Nhy+4,1)=-Nhy+1 bdc(Nhy+5,1)=-Nhy elseif(iresin.eq.3.or.iresin.eq.4)then !SER and CYS bdc(6,1)=-1 bdc(6,2)=Nhy-1 bdc(Nhy+5,1)=-Nhy+1 elseif(iresin.eq.5)then !VAL bdc(5,3)=2 bdc(6,1)=-1 bdc(6,2)=Nhy-2 bdc(6,3)=Nhy-1 bdc(6,3)=Nhy bdc(7,1)=-2 bdc(7,2)=Nhy bdc(7,3)=Nhy+1 bdc(7,4)=Nhy+2 bdc(Nhy+4,1)=-Nhy+2 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy bdc(Nhy+7,1)=-Nhy bdc(Nhy+8,1)=-Nhy-1 bdc(Nhy+9,1)=-Nhy-2 elseif(iresin.eq.6)then !THR bdc(5,2)=1 bdc(5,3)=2 bdc(6,1)=-1 bdc(6,2)=Nhy-2 bdc(7,1)=-2 bdc(7,2)=Nhy-2 bdc(7,3)=Nhy-1 bdc(7,4)=Nhy bdc(Nhy+4,1)=-Nhy+2 bdc(Nhy+5,1)=-Nhy+2 bdc(Nhy+6,1)=-Nhy+1 bdc(Nhy+7,1)=-Nhy elseif(iresin.eq.7)then !ILE bdc(5,3)=2 bdc(6,1)=-1 bdc(6,2)=2 bdc(6,3)=Nhy-2 bdc(6,4)=Nhy-1 bdc(7,1)=-2 bdc(7,2)=Nhy-1 bdc(7,3)=Nhy bdc(7,4)=Nhy+1 bdc(8,1)=-2 bdc(8,2)=Nhy+1 bdc(8,3)=Nhy+2 bdc(8,4)=Nhy+3 bdc(Nhy+4,1)=-Nhy+2 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy+1 bdc(Nhy+7,1)=-Nhy bdc(Nhy+8,1)=-Nhy-1 bdc(Nhy+9,1)=-Nhy-1 bdc(Nhy+10,1)=-Nhy-2 bdc(Nhy+11,1)=-Nhy-3 elseif(iresin.eq.8)then !PRO bdc(1,2)=Nhy-1 bdc(2,4)=Nhy-1 bdc(5,2)=1 bdc(5,3)=Nhy-3 bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=Nhy-2 bdc(6,4)=Nhy-1 bdc(7,1)=-1 bdc(7,2)=-Nhy+1 bdc(7,3)=Nhy-1 bdc(7,4)=Nhy bdc(Nhy+1,1)=-Nhy+1 bdc(Nhy+2,1)=-Nhy+3 bdc(Nhy+4,1)=-Nhy+2 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy+1 bdc(Nhy+7,1)=-Nhy elseif(iresin.eq.9)then !MET bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=Nhy-1 bdc(6,4)=Nhy bdc(7,1)=-1 bdc(7,2)=1 bdc(8,1)=-1 bdc(8,2)=Nhy-1 bdc(8,3)=Nhy bdc(8,4)=Nhy+1 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy bdc(Nhy+7,1)=-Nhy+1 bdc(Nhy+8,1)=-Nhy bdc(Nhy+9,1)=-Nhy-1 elseif(iresin.eq.10.or.iresin.eq.11)then !ASP and ASN bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=2 bdc(7,1)=-1 bdc(8,1)=-2 if(iresin.eq.11)then !ASN bdc(8,2)=Nhy-3 bdc(8,3)=Nhy-2 bdc(Nhy+5,1)=-Nhy+3 bdc(Nhy+6,1)=-Nhy+2 endif elseif(iresin.eq.12)then !LEU bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=2 bdc(6,4)=Nhy-1 bdc(7,1)=-1 bdc(7,2)=Nhy-1 bdc(7,3)=Nhy bdc(7,4)=Nhy+1 bdc(8,1)=-2 bdc(8,2)=Nhy+1 bdc(8,3)=Nhy+2 bdc(8,4)=Nhy+3 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy+1 bdc(Nhy+7,1)=-Nhy bdc(Nhy+8,1)=-Nhy-1 bdc(Nhy+9,1)=-Nhy-1 bdc(Nhy+10,1)=-Nhy-2 bdc(Nhy+11,1)=-Nhy-3 elseif(Iresin.eq.13)then !LYS bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=Nhy-1 bdc(6,4)=Nhy bdc(7,1)=-1 bdc(7,2)=1 bdc(7,3)=Nhy bdc(7,4)=Nhy+1 bdc(8,1)=-1 bdc(8,2)=1 bdc(8,3)=Nhy+1 bdc(8,4)=Nhy+2 bdc(9,1)=-1 bdc(9,2)=Nhy+2 bdc(9,3)=Nhy+3 bdc(9,4)=Nhy+4 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy bdc(Nhy+7,1)=-Nhy bdc(Nhy+8,1)=-Nhy-1 bdc(Nhy+9,1)=-Nhy-1 bdc(Nhy+10,1)=-Nhy-2 bdc(Nhy+11,1)=-Nhy-2 bdc(Nhy+12,1)=-Nhy-3 bdc(Nhy+13,1)=-Nhy-4 elseif(iresin.eq.14.or.iresin.eq.15)then !GLU & GLN bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=Nhy-1 bdc(6,4)=Nhy bdc(7,1)=-1 bdc(7,2)=1 bdc(7,3)=2 bdc(8,1)=-1 bdc(9,1)=-2 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy if(Iresin.eq.15)then bdc(9,2)=Nhy-2 bdc(9,3)=Nhy-1 bdc(Nhy+7,1)=-Nhy+2 bdc(Nhy+8,1)=-Nhy+1 endif elseif(Iresin.eq.16)then !ARG bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=Nhy-1 bdc(6,4)=Nhy bdc(7,1)=-1 bdc(7,2)=1 bdc(7,3)=Nhy bdc(7,4)=Nhy+1 bdc(8,1)=-1 bdc(8,2)=1 bdc(8,3)=Nhy+1 bdc(9,1)=-1 bdc(9,2)=1 bdc(9,3)=2 ! bdc(9,4)=1 bdc(10,1)=-1 bdc(10,2)=Nhy bdc(10,3)=Nhy+1 ! bdc(10,4)=-1 bdc(11,1)=-2 bdc(11,2)=Nhy+1 bdc(11,3)=Nhy+2 bdc(Nhy+5,1)=-Nhy+1 bdc(Nhy+6,1)=-Nhy bdc(Nhy+7,1)=-Nhy bdc(Nhy+8,1)=-Nhy-1 bdc(Nhy+9,1)=-Nhy-1 bdc(Nhy+10,1)=-Nhy bdc(Nhy+11,1)=-Nhy-1 bdc(Nhy+12,1)=-Nhy-1 bdc(Nhy+13,1)=-Nhy-2 elseif(Iresin.eq.17)then !HIS bdc(6,1)=-1 bdc(6,2)=1 bdc(6,3)=2 ! bdc(6,4)=2 bdc(7,1)=-1 bdc(7,2)=2 bdc(7,3)=Nhy-2 bdc(8,1)=-2 ! bdc(8,4)=-2 bdc(8,2)=2 bdc(8,3)=Nhy-2 bdc(9,1)=-2 bdc(9,2)=1 ! bdc(9,4)=1 bdc(9,3)=Nhy-2 bdc(Nhy+5,1)=-Nhy+2 bdc(Nhy+6,1)=-Nhy+2 bdc(Nhy+7,1)=-Nhy+2 elseif(Iresin.eq.18.or.Iresin.eq.19)then !PHE and TYR bdc(6,1)=-1 bdc(6,2)=1 ! bdc(6,4)=1 bdc(6,3)=2 bdc(7,1)=-1 ! bdc(7,4)=-1 bdc(7,2)=2 bdc(7,3)=Nhy-2 bdc(8,1)=-2 bdc(8,2)=2 ! bdc(8,4)=2 bdc(8,3)=Nhy-2 bdc(9,1)=-2 bdc(9,2)=2 ! bdc(9,4)=2 bdc(9,3)=Nhy-2 bdc(10,1)=-2 ! bdc(10,4)=-2 bdc(10,2)=1 bdc(10,3)=Nhy-2 bdc(11,1)=-2 bdc(11,2)=-1 ! bdc(11,4)=-2 bdc(11,3)=Nhy-2 bdc(Nhy+5,1)=-Nhy+2 bdc(Nhy+6,1)=-Nhy+2 bdc(Nhy+7,1)=-Nhy+2 bdc(Nhy+8,1)=-Nhy+2 bdc(Nhy+9,1)=-Nhy+2 if(Iresin.eq.19)then bdc(11,3)=1 bdc(11,4)=-2 bdc(12,1)=-1 bdc(12,2)=Nhy-3 bdc(Nhy+9,1)=-Nhy+3 endif elseif(Iresin.eq.20)then !TRP bdc(6,1)=-1 bdc(6,2)=1 ! bdc(6,4)=1 bdc(6,3)=2 bdc(7,1)=-1 bdc(7,2)=2 ! bdc(7,4)=-1 bdc(7,3)=Nhy-2 bdc(8,1)=-2 bdc(8,2)=2 ! bdc(8,4)=2 bdc(8,3)=3 bdc(9,1)=-2 bdc(9,2)=1 bdc(9,3)=Nhy-3 bdc(10,1)=-1 bdc(10,2)=-2 ! bdc(10,4)=-2 bdc(10,3)=2 bdc(11,1)=-3 bdc(11,2)=2 ! bdc(11,4)=2 bdc(11,3)=Nhy-4 bdc(12,1)=-2 bdc(12,2)=2 ! bdc(12,4)=2 bdc(12,3)=Nhy-4 bdc(13,1)=-2 bdc(13,2)=1 ! bdc(13,4)=-2 bdc(13,3)=Nhy-4 bdc(14,1)=-1 bdc(14,2)=-2 ! bdc(14,4)=-2 bdc(14,3)=Nhy-4 bdc(Nhy+5,1)=-Nhy+2 bdc(Nhy+6,1)=-Nhy+3 bdc(Nhy+7,1)=-Nhy+4 bdc(Nhy+8,1)=-Nhy+4 bdc(Nhy+9,1)=-Nhy+4 bdc(Nhy+10,1)=-Nhy+4 endif end subroutine Assbonda !!aasign the bond length value to 32*32 matrix subroutine bondp bd0(4,4)=1.335 !C-C bd0(23,23)=1.375 !CA-CA bd0(4,23)=1.49 !C-CA bd0(5,23)=1.49 !CC-CA bd0(30,31)=1.527 !CP1-CP2 bd0(31,31)=1.537 !CP2-CP2 bd0(31,32)=1.537 !CP2-CP3 bd0(6,6)=1.36 !CPH1-CPH1 bd0(8,23)=1.368 !CPT-CA bd0(8,8)=1.4 !CPT-CPT bd0(1,1)=1.5 !CT1-CT1 bd0(1,2)=1.538 !CT1-CT2 bd0(1,3)=1.538 !CT1-CT3 bd0(1,4)=1.49 !CT1-C bd0(1,5)=1.522 !CT1-CC bd0(2,2)=1.53 !CT2-CT2 bd0(2,3)=1.528 !CT2-CT3 bd0(2,4)=1.49 !CT2-C bd0(2,5)=1.522 !CT2-CC bd0(2,6)=1.5 !CT2-CPH1 bd0(2,14)=1.51 !CT2-CY bd0(2,23)=1.49 !CT2-CA bd0(3,3)=1.53 !CT3- bd0(3,4)=1.49 !CT3-C bd0(3,23)=1.49 !CT3-CA bd0(3,5)=1.522 !CT3-CC bd0(3,6)=1.5 !CT3-CPH1 bd0(14,23)=1.365 !CY-CA bd0(8,14)=1.44 !CPT-CY bd0(4,30)=1.49 bd0(5,30)=1.49 !!C finish bd0(18,23)=1.083 !HA-CA bd0(5,18)=1.1 !HA-CC bd0(18,31)=1.111 !HA-CP2 bd0(18,32)=1.111 !HA-CP3 bd0(1,18)=1.111 bd0(2,18)=1.111 bd0(3,18)=1.111 !HA-CT3 bd0(14,18)=1.08 !HA-CY bd0(19,31)=1.08 !HB-CP1 bd0(1,19)=1.08 bd0(2,19)=1.08 bd0(3,19)=1.08 !HB-CT3 bd0(23,25)=1.08 !HP-CA bd0(14,25)=1.08 !HP-CY bd0(6,22)=1.083 bd0(7,22)=1.09 !HR1-CPH2 bd0(19,30)=1.08 !!H finish bd0(4,9)=1.3 !N-C bd0(9,30)=1.434 !N-CP1 bd0(9,32)=1.455 !N-CP3 bd0(4,13)=1.365 !NC2-C bd0(2,13)=1.49 bd0(3,13)=1.49 !NC2-CT3 bd0(13,20)=1.0 !NC2-HC bd0(4,10)=1.345 bd0(1,10)=1.43 bd0(2,10)=1.43 bd0(3,10)=1.43 !NH1-CT3 bd0(10,17)=0.997 !NH1-H bd0(10,20)=0.98 !NH1-HC bd0(5,11)=1.36 !NH2-CC bd0(2,11)=1.455 bd0(3,11)=1.455 !NH2-CT3 bd0(1,12)=1.48 bd0(2,12)=1.48 bd0(3,12)=1.48 !NH3-CT3 bd0(12,20)=1.04 !NH3-HC bd0(6,15)=1.38 !NR2-CPH1 bd0(7,15)=1.32 bd0(16,23)=1.37 !NY-CA bd0(8,16)=1.375 !NY-CPT bd0(16,17)=0.976 !NY-H bd0(11,17)=1.0 bd0(15,17)=1.0 bd0(9,27)=1.455 !N finish bd0(4,26)=1.23 !O-C bd0(5,26)=1.23 bd0(23,27)=1.26 !OC-CA bd0(5,27)=1.26 !OC-CC bd0(4,27)=1.26 !OC-CC bd0(2,27)=1.33 bd0(3,27)=1.33 !OC-CT3 bd0(23,28)=1.411 !OH1-CA bd0(1,28)=1.42 bd0(2,28)=1.42 bd0(3,28)=1.42 !OH1-CT3 bd0(17,28)=0.96 !OH1-H bd0(2,29)=1.818 bd0(3,29)=1.816 !S-CT3 bd0(21,29)=1.325 !HS-S !1for PRO N_term bd0(30,33)=1.485 bd0(32,33)=1.502 bd0(20,33)=1.006 end subroutine bondp !!assign the value of bond angle subroutine angp ag0(23,23,23)=120.0 ag0(4,9,30)=117.0 ag0(4,30,31)=112.3 ag0(5,30,31)=112.3 ag0(30,31,31)=108.5 ag0(31,31,32)=108.5 ag0(4,9,32)=117.0 ag0(30,9,32)=114.2 ag0(6,15,7)=104.0 ag0(8,23,23)=118.0 ag0(8,8,23)=122.0 ag0(8,14,23)=107.4 ag0(8,16,23)=108.0 ag0(1,1,4)=108.0 ag0(1,1,5)=108.0 ag0(1,1,1)=111.0 ag0(1,1,2)=111.0 ag0(1,1,3)=108.5 ag0(1,2,1)=113.5 ag0(1,2,2)=113.5 ag0(1,2,3)=113.5 ag0(1,2,5)=108.0 ag0(1,2,6)=113.0 ag0(1,2,14)=114.0 ag0(1,2,23)=107.5 ag0(1,10,4)=120.0 ag0(2,23,23)=122.3 ag0(2,6,6)=130.0 ag0(2,1,3)=114.0 ag0(2,1,4)=108.0 ag0(2,1,5)=108.0 ag0(2,2,4)=108.0 ag0(2,2,5)=108.0 ag0(2,2,2)=113.5 ag0(2,2,3)=115.0 ag0(2,2,14)=114.0 ag0(2,14,23)=129.4 ag0(2,14,8)=124.0 ag0(2,13,4)=120.0 ag0(2,10,4)=120.0 ag0(2,29,3)=95.0 ag0(3,23,23)=122.3 ag0(3,6,6)=130.0 ag0(3,1,3)=114.0 ag0(3,1,4)=108.0 ag0(3,1,5)=108.0 ag0(3,1,23)=107.5 ag0(3,2,23)=107.5 ag0(3,2,6)=113.0 ag0(3,2,6)=114.0 ag0(3,13,4)=120.0 ag0(3,10,4)=109.6 ag0(8,8,14)=107.4 ag0(14,8,23)=130.0 ag0(4,10,17)=123.0 ag0(1,10,17)=117.0 ag0(2,10,17)=117.0 ag0(3,10,17)=117.0 ag0(5,10,17)=120.0 ag0(17,10,17)=120.0 ag0(17,16,23)=126.0 ag0(8,16,17)=126.0 ag0(17,28,23)=108.0 ag0(1,28,17)=106.0 ag0(2,28,17)=106.0 ag0(3,28,17)=106.0 ag0(18,23,23)=120.0 ag0(8,23,18)=122.0 ag0(14,23,18)=125.0 ag0(18,31,30)=110.0 ag0(18,31,31)=110.0 ag0(18,31,32)=110.0 ag0(18,31,18)=109.0 ag0(18,32,31)=110.0 ag0(18,32,18)=109.0 ag0(4,1,18)=109.5 ag0(1,1,18)=110.0 ag0(2,1,18)=110.0 ag0(3,1,18)=110.0 ag0(18,1,18)=109.0 ag0(4,2,18)=109.5 ag0(5,2,18)=109.5 ag0(6,2,18)=109.5 ag0(18,2,23)=107.5 ag0(1,2,18)=110.0 ag0(2,2,18)=110.0 ag0(3,2,18)=110.0 ag0(14,2,18)=109.5 ag0(18,2,18)=109.5 ag0(4,3,18)=109.5 ag0(5,3,18)=109.5 ag0(6,3,18)=109.5 ag0(18,3,23)=107.5 ag0(1,3,18)=110.0 ag0(2,3,18)=110.0 ag0(3,3,18)=110.0 ag0(18,3,18)=109.4 ag0(18,16,23)=126.4 ag0(8,16,18)=126.4 ag0(4,30,19)=112.0 ag0(5,30,19)=112.0 ag0(4,1,19)=109.5 ag0(5,1,19)=109.5 ag0(1,1,19)=111.0 ag0(2,1,19)=111.0 ag0(3,1,19)=111.0 ag0(4,2,19)=109.5 ag0(5,2,19)=109.5 ag0(19,2,19)=115.0 ag0(4,3,19)=109.5 ag0(2,13,20)=120.0 ag0(3,13,20)=120.0 ag0(4,13,20)=120.0 ag0(20,13,20)=120.0 ag0(2,11,20)=111.0 ag0(3,11,20)=111.0 ag0(20,11,20)=106.5 ag0(1,12,20)=109.5 ag0(2,12,20)=109.5 ag0(3,12,20)=109.5 ag0(20,12,20)=109.5 ag0(23,12,25)=120.0 ag0(8,12,25)=122.0 ag0(14,12,25)=125.0 ag0(23,14,25)=126.4 ag0(8,14,25)=126.4 ag0(6,6,22)=130.0 ag0(6,6,24)=130.0 ag0(5,11,17)=120.0 ag0(7,15,17)=126.0 ag0(15,7,22)=125.0 ag0(23,23,25)=120.0 ag0(2,29,21)=95.0 ag0(3,29,21)=95.0 !!H finish ag0(9,4,30)=112.5 ag0(1,4,9)=112.5 ag0(2,4,9)=112.5 ag0(3,4,9)=112.5 ag0(4,30,9)=108.2 ag0(5,30,9)=108.2 ag0(9,30,31)=110.8 ag0(9,30,19)=112.0 ag0(9,32,18)=108.0 ag0(9,32,31)=110.5 ag0(13,4,13)=120.0 ag0(2,2,13)=107.5 ag0(13,2,18)=107.5 ag0(13,3,18)=107.5 ag0(10,4,30)=116.5 ag0(1,4,10)=116.5 ag0(2,4,10)=116.5 ag0(3,4,10)=116.5 ag0(1,1,10)=113.5 ag0(2,1,10)=113.5 ag0(3,1,10)=113.5 ag0(4,1,10)=107.0 ag0(5,1,10)=107.0 ag0(10,1,19)=108.0 ag0(2,2,10)=113.5 ag0(4,2,10)=107.0 ag0(5,2,10)=107.0 ag0(10,2,18)=109.5 ag0(10,2,19)=108.0 ag0(10,3,18)=109.5 ag0(11,5,30)=112.5 ag0(1,5,11)=116.5 ag0(2,5,11)=116.5 ag0(3,5,11)=116.5 ag0(11,5,18)=111.0 ag0(11,2,18)=109.5 ag0(11,2,19)=109.5 ag0(2,2,11)=110.0 ag0(11,3,18)=109.5 ag0(1,1,12)=110.0 ag0(2,1,12)=110.0 ag0(3,1,12)=110.0 ag0(4,1,12)=110.0 ag0(5,1,12)=110.0 ag0(12,1,19)=107.5 ag0(4,2,12)=110.0 ag0(5,2,12)=110.0 ag0(2,2,12)=110.0 ag0(12,2,18)=107.5 ag0(12,2,19)=107.5 ag0(12,3,18)=107.5 ag0(6,6,15)=110.0 ag0(2,6,15)=120.0 ag0(15,6,24)=120.0 ag0(15,6,22)=125.0 ag0(14,23,16)=110.0 ag0(16,23,18)=125.0 ag0(16,23,25)=125.0 ag0(16,8,23)=130.6 ag0(8,8,16)=107.4 !!N finish ag0(26,4,30)=118.0 ag0(1,4,26)=121.0 ag0(2,4,26)=121.0 ag0(3,4,26)=121.0 ag0(17,4,26)=121.0 ag0(9,4,26)=122.5 ag0(10,4,26)=122.5 ag0(26,5,30)=118.0 ag0(1,5,26)=121.0 ag0(2,5,26)=121.0 ag0(3,5,26)=121.0 ag0(18,5,26)=122.0 ag0(11,5,26)=122.5 ag0(23,23,27)=120.0 ag0(27,5,30)=118.0 ag0(1,5,27)=118.0 ag0(2,5,27)=118.0 ag0(3,5,27)=118.0 ag0(27,5,27)=124.0 ag0(3,2,27)=122.0 ag0(18,2,27)=118.3 ag0(19,2,27)=118.3 ag0(23,23,28)=120.0 ag0(1,1,28)=110.1 ag0(3,1,28)=110.1 ag0(18,1,28)=108.89 ag0(1,2,28)=110.1 ag0(2,2,28)=110.1 ag0(3,2,28)=110.1 ag0(26,5,27)=124.0 ag0(18,2,28)=108.89 ag0(18,3,28)=108.89 !!O finish ag0(1,2,29)=112.5 ag0(2,2,29)=114.5 ag0(3,2,29)=114.5 ag0(18,2,29)=111.3 ag0(18,3,29)=111.3 !!for PRO N_term ag0(30,33,32)=111.0 ag0(17,33,31)=109.5 ag0(17,33,33)=109.5 ag0(17,33,17)=107.5 ag0(4,30,33)=106.0 ag0(5,30,33)=106.0 ag0(32,30,33)=108.5 ag0(19,30,33)=107.5 ag0(30,32,33)=108.5 ag0(18,32,33)=109.15 ag0(20,33,32)=109.5 end subroutine angp end program main