! Reads coordinates in a pdb format subroutine rdpdb(ierror) ! IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" character List_res*200,List_lig*1000,Line*120,res_file*60, . digit*10,List_pres*6,at_symb*2,lc_list*26,up_list*26, . List_nres*6,old_rnum*4 logical found_res,iex,first,dna,found_at,protein,end_of_list .,res_pointer, ligand integer charge,atom_idx data digit/'0123456789'/ data List_pres/'LYSARG'/ data List_nres/'ASPGLU'/ data lc_list/'abcdefghijklmnopqrstuvwxyz'/ data up_list/'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ !-------------------------------------------------------------------- ! Check if there is a pdb file inquire( file=fname(5), exist=iex) if(.not.iex) then C-RDC write(iscr,*) 'rdpdb: pdb file not found ... EXIT !' ierror = 1 return else ! Open pdb file open(ipdb,file=fname(5),status='old') endif ! Read residue labels, open file with residue/ligand names c i=iargc() res_file = "residue.dat" c if(i.gt.0) call getarg(1,res_file) inquire( file=res_file, exist=iex) if(.not.iex) then C-RDC write(iscr,*) 'rdpdb: residue file not found ... EXIT !' ierror = 1 return else ! Open file with residue labels open(50,file=res_file,status='old') ! The file with residue labels is divided in two parts: amino acids ! residues or DNA nucleotides are listes in the upper part which ends ! with the separator '---'. The lower part is reserved for ligands. ! This division is used to recognize if a molecule is of DNA, PROTEIN ! or LIGAND type and this information is used to set up pointers for ! subsystems in the following way: ! for 1 - DNA the pointer is set on the O of the poshphate group ! 2 - PROTEIN the pointer is set on the backbone N of the residue ! 3 - LIGAND the pointer is set on the first atom listed in pdb file ! and all other atoms for which the keyword RES is in the input line ! n_res = -1 i=1 end_of_list = .false. do while (.not.end_of_list) read(50,'(a3)',End=3) List_res(i:i+2) if(List_res(i:i+2).eq.'---') end_of_list = .true. C-RDC ! write(iscr,'(a,a3)') 'DNA/PROT: ',List_res(i:i+2) n_res = n_res + 1 i = i + 3 enddo ! Read 3-letter ligand labels n_lig = 0 i=1 List_lig(1:3) = ' ' do while (.true.) read(50,'(a3)',End=4) List_lig(i:i+2) C-RDC ! write(iscr,'(a,a3)') 'LIGAND: ',List_lig(i:i+2) n_lig = n_lig + 1 i = i + 3 enddo 4 close(50) endif ! atom info: ! atom = 1 :5 ! res = 7 :10 ! chn = 11:11 ! rnum = 12:15 l = 0 k = 1 ires = 0 mol = 1 I = 0 charge = 0 nres = 0 old_rnum = ' ' do while (.True.) l = l + 1 line(1:6) = '@@@@@@' read(ipdb,'(a120)',End=1) Line if(line(1:6).eq.'ATOM'.or.line(1:6).eq.'HETATM') then IF(I.EQ.MAXATM)THEN IERROR = 1 C-RDC WRITE(IOUT,120) 120 FORMAT(/' MAXIMUM ALLOWED NUMBER OF ATOMS REACHED', . /' -- INCREASE MAXATM PARAMETER IN divcon.dim', . ' AND RECOMPILE') RETURN ENDIF I = I + 1 atom_inf(I)(1:15) = '##### #### ####' C-RDC ! write(iout,'(i6,3x,a80)') l,Line C-RDC ! Slice the one line of the pdb file istart = 7 ! skip ATOM or HETATM words call rdword(line,istart,istop) call rdinum(line,istart,atom_idx,ierror) ! read atom index IF(ierror.ne.0) then C-RDC write(iscr,*)'rdpdb: reading atom index in line ',l ierror = 1 return endif ISTART = ISTOP + 1 call rdword(line,istart,istop) ! read atom label ! If atom symbol starts with a digit pass the digit to the end of the label if(index(digit,line(ISTART:ISTART)).ne.0) then atom_inf(I)(1:5) = line(ISTART+1:ISTOP) LAST =ISTOP-ISTART + 1 atom_inf(I)(LAST:LAST) = line(ISTART:ISTART) line(ISTART:ISTOP) = atom_inf(I)(1:5) else atom_inf(I)(1:5) = line(ISTART:ISTOP) endif ! Extract the atom symbol from atom label at_symb = ' ' if(index(lc_list,atom_inf(I)(2:2)).ne.0) then ! Atom symbol made of 2 letter at_symb = atom_inf(I)(1:2) else ! Atom symbol made of 1 letter at_symb = atom_inf(I)(1:1) endif CALL UPCASE1(at_symb,2) j = 0 found_at = .false. DO While (.not.found_at.and.j.le.82) IF(at_symb.EQ.SYMBOL(J))THEN IATNUM(I) = J found_at = .true. ENDIF j = j + 1 ENDDO IF(.not.found_at) then IERROR = 1 C-RDC WRITE(IOUT,'(/" UNRECOGNIZED ELEMENTAL SYMBOL: ",A2)') C-RDC . at_symb RETURN ENDIF ISTART = ISTOP + 1 call rdword(line,istart,istop) ! read residue label atom_inf(I)(7:10) = line(istart:istop) found_res = .false. dna = .false. protein = .false. if( istart.eq.istop ) then ! This may be a 1 letter DNA residue if( index(List_res,atom_inf(I)(7:7)).eq.0 ) Then ! DNA residue is not in the list C-RDC Write(iscr,'(a,a1,a,i4)') C-RDC . 'rdpdb: Unmatched 1 letter DNA residue (',atom_inf(I)(7:7), C-RDC . ') in line ',l ierror = 1 return else dna = .true. ! 1 letter DNA residue found in the list found_res = .true. ! Move letter from 'third' to 'first' position as required by AMBER line(istart-2:istart-2) = line(istart:istart) line(istart:istart) = ' ' endif elseif(istop-istart.eq.1) then ! This may be a 2 letter DNA residue if( index(List_res,atom_inf(I)(7:8)).eq.0 ) Then ! DNA residue is not in the list C-RDC Write(iscr,'(a,a2,a,i4)') C-RDC . 'rdpdb: Unmatched 2 letter DNA residue (', C-RDC . atom_inf(I)(7:8),') in line ',l ierror = 1 return else dna = .true. ! 2 letter DNA residue found in the list found_res = .true. ! Move letter from 'second' to 'first' and 'third' to 'second' position as required by AMBER line(istart-1:istart-1) = line(istart:istart) line(istart:istart) = line(istop:istop) line(istart:istart) = ' ' endif elseif( index(List_res,atom_inf(I)(7:9)).ne.0 ) then ! 3 letter amino acid residue found in the list found_res = .true. protein = .true. elseif( index(List_lig,atom_inf(I)(7:9)).ne.0 ) then ! 3 letter ligand found in the list found_res = .true. ligand = .true. else ! This residue label is not in the list of standard residue labels C-RDC Write(iscr,'(a,a4,a,i4)') C-RDC . 'rdpdb: Unmatched 3 letter residue/ligand (',atom_inf(I)(7:10), C-RDC . ') in line ',l ierror = 1 return endif ISTART = ISTOP + 1 ! Look for the residue tag/number call rdword(line,istart,istop) j = 12 if(index(up_list,line(istart:istart)).ne.0) then ! The residue number has a tag to distinguish between several proteins in the ! same pdb file atom_inf(I)(12:12) = line(istart:istart) ISTART = ISTART + 1 ! Look for the residue number call rdword(line,istart,istop) j = 13 endif atom_inf(I)(j:15) = line(istart:istop) ISTART = ISTOP + 1 ! Read cartesian coordinates call rdword(line,istart,istop) call getnum(line,istart,istop,xyz(1,I),ierror) ! read X IF(ierror.ne.0) goto 2 ISTART = ISTOP + 1 call rdword(line,istart,istop) call getnum(line,istart,istop,xyz(2,I),ierror) ! read Y IF(ierror.ne.0) goto 2 ISTART = ISTOP + 1 call rdword(line,istart,istop) call getnum(line,istart,istop,xyz(3,I),ierror) ! read Z IF(ierror.ne.0) goto 2 IF(RESDUE)THEN ! setup pointer on DNA nucleotide res_pointer = .false. if(dna.and.(atom_inf(I)(1:2).eq.'P ' . .or.atom_inf(I)(1:4).eq."O5' " . .or.atom_inf(I)(1:4).eq."O5* " . .or.index(line,'RES').ne.0)) then res_pointer = .true. if(old_rnum.ne.atom_inf(I)(12:15)) then ! sum up charges on DNA nucleotides old_rnum = atom_inf(I)(12:15) charge = charge - 1 endif elseif(protein.and.(atom_inf(I)(1:2).eq.'N ' . .or.index(line,'RES').ne.0) ) then ! setup pointer on amino acid residue res_pointer = .true. if(old_rnum.ne.atom_inf(I)(12:15)) then ! sum up charges on charged amino acid residues old_rnum = atom_inf(I)(12:15) if(index(List_nres,atom_inf(I)(7:9)).ne.0) then ! that's a negatively charged residue charge = charge - 1 C-RDC ! write(iscr,'(a15,a)') atom_inf(I)(1:15),' -RES' elseif(index(List_pres,atom_inf(I)(7:9)).ne.0 ) then ! that's a positively charged residue charge = charge + 1 C-RDC ! write(iscr,'(a15,a)') atom_inf(I)(1:15),' +RES' endif endif elseif(old_rnum.ne.atom_inf(I)(12:15) . .or.index(line,'RES').ne.0) then ! setup pointer on ligand for first atom and those with RES keyword res_pointer = .true. if(old_rnum.ne.atom_inf(I)(12:15)) . old_rnum = atom_inf(I)(12:15) endif if(res_pointer) then NRES = NRES + 1 IF(NRES.EQ.MAXRES)THEN IERROR = 1 C-RDC WRITE(IOUT, C-RDC . '(/" MAXIMUM ALLOWED NUMBER OF RESIDUES REACHED", C-RDC . /" -- INCREASE MAXRES PARAMETER IN divcon.dim", C-RDC . " AND RECOMPILE")') RETURN ENDIF IRPNT(NRES) = I C-RDC ! write(iout,'(i6,5x,a15,a)') atom_idx,atom_inf(I)(1:15) endif ENDIF elseif(line(1:3).eq.'TER') then ! Molecule finished k = 1 C-RDC write(iscr,*) 'Molecule :',mol,', residues: ',ires mol = mol + 1 endif C-RDC ! write(iout,'(i6,5x,a15)') i,atom_inf(i) enddo 1 close(ipdb) natoms = I C C END OF COORDINATE DATA. DO SOME CHECKING. C IF(RESDUE)THEN IF(NRES.EQ.0)THEN IERROR = 1 C-RDC WRITE(IOUT, C-RDC . '(/" ''RESIDUE'' KEYWORD SPECIFIED, BUT NO", C-RDC . " RESIDUES WERE DEFINED")') RETURN ENDIF C C ADD POINTER FOR NRES+1 FOR CONVENIENCE IN LOOPING. C IRPNT(NRES+1) = NATOMS + 1 ENDIF C-RDC write(iscr,'(2(a,i6),a,i3)') 'rdpdb: nres = ',nres, C-RDC . ', natoms = ',natoms,', total charge=',charge ! stop return 2 continue C-RDC 2 write(iscr,'(a)') C-RDC . 'rdpdb: Error in reading coordinates in pdb file' ierror = 1 return 3 continue C-RDC 3 write(iscr,'(a)') 'rdpdb: Unexpected EOF in residue.dat' ierror = 1 return end