program mardi2amber cccccccccccccccccccccccccccccccc c Current version 2.0 - Last change February 6th, 2001 c c written by He Liu at UCSF c Jul 2 1992 c c re-written by Marco Tonelli at UCSF c Jan 18 2001 c c write a NMR distance constraint file for Amber c (namelist format only) c c r1 is set to 0.0 c r2, r3 are the lower or upper bands (average distance -/+ SD). c r4 = r3 + rs (rs is the width of the parabolic wall) c c the force constant k1 (= k2) are the same for all constraints. c c INPUT FILES : c constraint file with format: c HA21 2 HA22 2 3.2 4.4 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c New Feb 2001 : c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c if the format of the constraint file is different, it is c easy to adapt the program to the different format, all you c need to do is change iat1,ire2,iat2,ire2 that specify the c column positions for atom names and residue numbers, read c below for more detail c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c c lower and upper bound can be any column after atom names and c residue numbers (their position must be specified) c c a pdb file which gives the absolute numbers or the proton c names for pseudo atoms. c the atom number following the ATOM keyword is used (so the c user needs to make sure that it is correct) c c To compile this program type the following at the prompt (you need c to have a fortran compiler): c c f77 -o mardi2amber mardi2amber.f c cccccccccccccccccccccccccccccccccc parameter(mat=1000,MW=40) parameter(len=200) character*79 fname character*(len) line character*5 atnm real rb(4) real fc,rs integer Position integer ia1(3),ia2(3) integer ic1,ic2 integer ib(MW),ie(MW) integer ntf(40) c pdb file character*5 patm(mat) integer Inmb(mat),Ires(mat) c constraint file character*5 at1,at2 integer re1,re2 real dlow,dup i0=ichar('0') iw=MW c c unit 1 = constraint file c unit 3 = pdb file c unit 4 = output file c icnstr=1 ipdb=3 iout=4 c c ic1 column position with lower bound c ic2 column position with upper bound c fc force constant c rs well width c c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c NEW February 6 2001 c CONSTRAINT FILE, ATOM NAMES FORMAT : c iat1,ire1,iat2,ire2 specify the column position for c (iat1) name and (ire1) residue # for atom 1 c (iat2) name and (ire2) residue # for atom 2 c c for example if the constraint file has a format c like this : c c HA1 1 HB1 1 3.6 4.5 c HA1 2 HB1 2 2.6 5.3 c HA1 3 HB1 3 3.4 4.2 c ^ ^ ^ ^ c (column # : 1 2 3 4) c c set values to : iat1=1 , ire1=2 , iat2=3 , ire2=4 c and compile the program c iat1 = 1 ire1 = 2 iat2 = 3 ire2 = 4 c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c c c 10 write (6,'(/"Enter input constraint file: ",$)') read (5,'(a)') fname open(icnstr,file=fname,status='old',err=10) 13 write (6,'(/"Enter number for position of LOWER bound column: ",$)') read (5,*,err=13) ic1 16 write (6,'(/"Enter number for position of UPPER bound column: ",$)') read (5,*,err=16) ic2 20 write (6,'(/"Enter pdb filename: ",$)') read (5,'(a)') fname open(ipdb,file=fname,status='old',err=20) write (6,'(/"Enter output file: ",$)') read (5,'(a)') fname open(iout,file=fname,status='unknown') 23 write (6,'(/"Enter force constant: ",$)') read (5,*,err=23) fc 26 write (6,'(/"Enter width of the parabolic wall: ",$)') read (5,*,err=26) rs c$$$ write (6,'(/a,$)') 'Select format: (n) for namelist, (i) for interface: ' c$$$ read (5,'(a1)') fmt c$$$ lstr=8 c$$$ if(fmt.eq.'n'.or.fmt.eq.'N') lstr=4 c c read pdb file : c - read only H atoms c - allow for two different formats (the difference is in the position c of the atom name column that is shifted by one) c np=0 50 read (ipdb,'(a79)',end=60) fname if(fname(1:4).eq.'ATOM'.and.(fname(13:13).eq.'H'.or. 1 fname(13:14).eq.'1H'.or.fname(13:14).eq.'2H'.or. 2 fname(13:14).eq.'3H'.or.fname(14:15).eq.'1H'.or. 3 fname(14:15).eq.'3H'.or.fname(14:15).eq.'2H'.or. 4 fname(14:14).eq.'H')) then np=np+1 if(fname(13:13).ne.' ') then read(fname,'(5x,i6,1x,a4,7x,i3)') Inmb(np),patm(np),Ires(np) else read(fname,'(5x,i6,2x,a4,6x,i3)') Inmb(np),patm(np),Ires(np) endif endif goto 50 60 print '(/"> Number of protons found in pdb file: ",i5)', np c c read constraint file c idst=0 100 read (icnstr,'(a)',end=110) line if (line(1:1).eq.'#') goto 100 call words(line,len,il,ib,ie,iw) read (line(ib(iat1):ie(iat1)),'(a)') at1 read (line(ib(ire1):ie(ire1)),*,err=100) re1 read (line(ib(iat2):ie(iat2)),'(a)') at2 read (line(ib(ire2):ie(ire2)),*,err=100) re2 read (line(ib(ic1):ie(ic1)),*,err=100) dlow read (line(ib(ic2):ie(ic2)),*,err=100) dup c c check that lower bounds is actually smaller than upper bound c otherwise skip line c if (dlow.gt.dup) then print '(/"> WARNING: Lower bound is BIGGER than Upper bound")' print '( "> Distance >>",a4,i3,x,a4,i3,"<< skipped.")', at1,re1,at2,re2 goto 100 endif c c set nps# : c 1 for regular H atoms c 2 for pseudoatoms of type Q or R c 3 for pseudoatoms of type M c nps1=1 nps2=1 if (at1(1:1).eq.'M') then nps1=3 elseif(at1(1:1).eq.'Q'.or.at1(1:1).eq.'R') then nps1=2 endif if (at2(1:1).eq.'M') then nps2=3 elseif(at2(1:1).eq.'Q'.or.at2(1:1).eq.'R') then nps2=2 endif c c search ATOM1 in list of atoms from pdb file c i1=0 if (nps1.eq.1) then ii = Position(at1,re1,patm,Ires,np) if(ii.gt.0) then ia1(1)=ii i1=i1+1 endif else lat = INDEX(at1," ")-1 do n=1,nps1 atnm = "H"//at1(2:lat)//char(i0+n) ii = Position(atnm,re1,patm,Ires,np) if (ii.eq.0) then atnm = char(i0+n)//"H"//at1(2:) ii = Position(atnm,re1,patm,Ires,np) endif if (ii.gt.0) then ia1(n)=ii i1=i1+1 endif enddo endif c c check if atom(s) 1 was(were) found c if (i1.ne.nps1) then print '(/"> WARNING: Atom >>",a4,i3,"<< not found in pdb file,")', at1,re1 print '( "> Distance >>",a4,i3,x,a4,i3,"<< skipped.")', at1,re1,at2,re2 goto 100 endif c c search ATOM2 in list of atoms from pdb file c i2=0 if (nps2.eq.1) then ii = Position(at2,re2,patm,Ires,np) if (ii.gt.0) then ia2(1)=ii i2=i2+1 endif else lat = INDEX(at2," ")-1 do n=1,nps2 atnm = "H"//at2(2:lat)//char(i0+n) ii = Position(atnm,re2,patm,Ires,np) if (ii.eq.0) then atnm = char(i0+n)//"H"//at2(2:) ii = Position(atnm,re2,patm,Ires,np) endif if (ii.gt.0) then ia2(n)=ii i2=i2+1 endif enddo endif c c check if atom(s) 2 was(were) found c if (i2.ne.nps2) then print '(/"> WARNING: Atom >>",a4,i3,"<< not found in pdb file.")', at2,re2 print '( "> Distance >>",a4,i3,x,a4,i3,"<< skipped.")', at1,re1,at2,re2 goto 100 endif c$$$ print *, (patm(ia1(j)),Ires(ia1(j)),j=1,i1) c$$$ print *, (patm(ia2(j)),Ires(ia2(j)),j=1,i2) c c c WRITE AMBER OUTPUT FILE c c c set bounds for amber c c uncomment these lines if you wish to set r1 = r2 - well width c$$$ rb(1)=dlow-rs c$$$ if(rb(1).lt.0.) rb(1)=0. rb(1)=0.0 rb(2)=dlow rb(3)=dup rb(4)=rb(3)+rs c c &rst iat= 151, 1254, c r1= 0.0, r2= 0.0, r3= 3.50, r4= 5.00, c rk2= 10.00, rk3= 10.00 &end c idst=idst+1 write (iout,999) idst,at1,re1,at2,re2 if (nps1.eq.1 .AND. nps2.eq.1) then write (iout,1000) Inmb(ia1(1)),Inmb(ia2(1)) elseif (nps1.gt.1 .AND. nps2.eq.1) then write (iout,1001) (Inmb(ia1(i)),i=1,i1),Inmb(ia2(1)) elseif (nps1.eq.1 .AND. nps2.gt.1) then write (iout,1002) Inmb(ia1(1)),(Inmb(ia2(i)),i=1,i2) elseif (nps1.gt.1 .AND. nps2.gt.1) then write (iout,1003) (Inmb(ia1(i)),i=1,i1),(Inmb(ia2(i)),i=1,i2) endif write (iout,1005) (i,rb(i),i=1,4),fc,fc 999 format('# Distance constraint number',i5,' : ',2(a4,i3,x)) 1000 format(' &rst iat=',i4,',',i4',') 1001 format(' &rst iat= -1,',i4,',',/t7'igr1 =',(i4,',')) 1002 format(' &rst iat=',i4,', -1,',/t7'igr2 =',(i4,',')) 1003 format(' &rst iat= -1, -1,' 1 ,/t7' igr1 =',(i4,','),' igr2 =',(i4,',')) 1005 format(t7,4('r',i1,'=',f6.2,', '), 1 /t7,'rk2=',f7.2,', rk3=',f7.2,' &end') goto 100 110 continue print '(/"> Number of constraints written to output file:",i5)', idst stop end c c INTEGER function Position : c returns an integer that corresponds to the position of the atom c in the atom list from pdb file or 0 if the atom is not found c INTEGER function Position (at,rs,atom,residue,itot) character*(*) at,atom(*) integer rs,residue(*),itot n=0 do i=1,itot if (at.eq.atom(i) .AND. rs.eq.residue(i)) then n=i goto 100 endif enddo 100 Position=n return end subroutine words(rec,len,m,ib,ie,k) c Anil Kumar December 3, 1992 c modified Marco Tonelli February 5, 2001 character*(*) rec integer ib(*),ie(*) logical w w=.true. do i=1,k ib(i)=0 ie(i)=0 enddo k=0 do i=1,len if((rec(i:i).ne.' ').and. 1 (ichar(rec(i:i)).ne.9))then if(w)then k=k+1 ib(k)=i endif ie(k)=i w=.false. else w=.true. endif enddo m=ie(k) return end