Commit d52100bb authored by Steffen Möller's avatar Steffen Möller
Browse files

New upstream version 20190822+dfsg

parent d848d030

TMalign.cpp

0 → 100644
+5137 −0

File added.

Preview size limit exceeded, changes collapsed.

TMalign.f

+94 −35
Original line number Diff line number Diff line
@@ -77,6 +77,8 @@ 
*     2017/07/08: Added one iteration in initial4 to avoid asymmetric alignment

*     2019/07/08: Enable TM-align to support both PDB and mmCIF formats, and

*                 fixed a bug on file output

*     2019/08/18: Fixed multiple bugs associated with mmCIF formats

*     2019/08/22: added output scripts for pymol, C++ version was included.

**************************************************************************

      

c        1         2         3         4         5         6         7 !
@@ -102,7 +104,7 @@ c 3456789012345678901234567890123456789012345678901234567890123456789012345678 
      character*500 fnam,pdb(100),outname,falign,fmatrix

      character*3 aa(-1:20),aanam,ss(2,nmax)

      character*500 s,du,dum1,dum2

      character*504 outnameall_tmp,outnameall

      character*504 xxx,xxx_p,xxx_pdb

      character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2)

      character*8 version

      character*5000 s1         !maximum length of protein is 5000
@@ -113,11 +115,13 @@ ccc mmCIF 
      character*20 Aatom(2,90000)

      character Agroup(2,90000)*6

      character Ares(2,90000)*3,Aalt(2,90000)

      character Ains(2,9000),Ach(2,90000),Aent(2,90000)

      character Ains(2,90000),Aent(2,90000)

      character Cins(2,nmax),Cch(2)

      character*2 Ach(2,90000),chid(2,90000)

      

      integer Aatomi(2,90000),Aresi(2,90000),n_cut(2)

      dimension iform(10),xx(2,90000),yy(2,90000),zz(2,90000),nL(2)

      dimension iform(10),xx(2,90000),yy(2,90000),zz(2,90000)

      dimension nL(2),nLCA(2)

ccc

      

      dimension m1(nmax),m2(nmax),m12(2,nmax)
@@ -179,20 +183,25 @@ ccc 
     &        '''TM.sup_all'' and ''TM.sup_atm'':'

         write(*,*)'   >TMalign chain_1.pdb chain_2.pdb -o TM.sup'

         write(*,*)'      To view superimposed C-alpha traces of',

     &        ' aligned regions by rasmol:'

     &        ' aligned regions by rasmol or pymol:'

         write(*,*)'        >rasmol -script TM.sup'

         write(*,*)'        >pymol -d @TM.sup.pml'

         write(*,*)'      To view superimposed C-alpha traces of',

     &        ' all regions:'

         write(*,*)'        >rasmol -script TM.sup_all'

         write(*,*)'        >pymol -d @TM.sup_all.pml'

         write(*,*)'      To view superimposed full-atom structures of',

     &        ' aligned regions:'

         write(*,*)'        >rasmol -script TM.sup_atm'

         write(*,*)'        >pymol -d @TM.sup_atm.pml'

         write(*,*)'      To view superimposed full-atom structures of',

     &        ' all regions:'

         write(*,*)'        >rasmol -script TM.sup_all_atm'

         write(*,*)'        >pymol -d @TM.sup_all_atm.pml'

         write(*,*)'      To view superimposed full-atom structures of',

     &        ' all regions with ligands:'

         write(*,*)'        >rasmol -script TM.sup_all_atm_lig'

         write(*,*)'        >pymol -d @TM.sup_all_atm_lig.pml'

         write(*,*)

         write(*,*)'4. There are two TM-scores reported. You ',

     &        'should use the one normalized by'
@@ -214,7 +223,7 @@ ccc 
         goto 9999

      endif

      

      version='20190708'

      version='20190822'

      if(fnam.eq.'-v')then

         write(*,*)'TM-align Version ',version

         goto 9999
@@ -333,8 +342,8 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^

                     if(mk.eq.1)then !mk=1 ------------>

                        i=i+1

                        read(s,'(a6,I5,a6,A3,A2,i4,A1,a3,3F8.3)')

     &                       du,itmp,du,aanam,du,mm(ic,i),

                        read(s,'(a6,I5,a6,A3,A1,A1,i4,A1,a3,3F8.3)')

     &                       du,itmp,du,aanam,du,chid(ic,i),mm(ic,i),

     &                       Cins(ic,i),du,

     &                       xa(1,i,ic-1),xa(2,i,ic-1),xa(3,i,ic-1)

***
@@ -363,8 +372,10 @@ c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^ 
                  if(du.eq.'_atom_site.label_alt_id')i_alt=in !'.',A,B

                  if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU

                  if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B

                  if(du.eq.'_atom_site.auth_asym_id')i_ch=in !A,B, using later one

                  if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b

                  if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3

                  if(du.eq.'_atom_site.auth_seq_id')i_resi=in !1,2,3, identical to PDB res

                  if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?

                  if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234

                  if(du.eq.'_atom_site.Cartn_y')i_y=in !y, 1.234
@@ -410,6 +421,7 @@ c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^ 
                        ss(ic,i)=ctmp(i_res) !residue name, 'GLY', for seq_ID

***   

                        i8=mm(ic,i)

                        chid(ic,i)=ch(i)

                        nres1(ic,i8,ichar(du1))=

     &                       nres1(ic,i8,ichar(du1))+1 !nres1 only for check altLoc

***
@@ -802,6 +814,7 @@ c************* Step-1: read all-atom structures ****************** 
c*****************************************************************

      do 1002 ic=1,2            !ic=1,2 for file1 and file2

         nL(ic)=0               !number of lines in PDB file

         nLCA(ic)=0             !number of lines with CA in PDB file

         n_cut(ic)=0            !end of first chain, decide to show aaa_all_atm_lig

         open(unit=10,file=pdb(ic),status='old')

         if(iform(ic).eq.1)then !file in PDB format %%%%%%%%%%------->
@@ -845,7 +858,8 @@ c^^^^^^^^^altLoc checked (mk.ne.1) will be skipped) ^^^^^^^^^ 
                  

                  if(mk.eq.1)then !mk=1--------->

                     nL(ic)=nL(ic)+1

                     read(s,8999)Agroup(ic,nL(ic)),Aatomi(ic,nL(ic)),

                     read(s,'(a6,I5,a1,A4,a1,A3,a1,A1,I4,A1,a3,3F8.3)')

     &                    Agroup(ic,nL(ic)),Aatomi(ic,nL(ic)),

     &                    du,Aatom(ic,nL(ic)),Aalt(ic,nL(ic)),

     &                    Ares(ic,nL(ic)),du,Ach(ic,nL(ic)),

     &                    Aresi(ic,nL(ic)),Ains(ic,nL(ic)),du,
@@ -883,8 +897,10 @@ c read(Aatom(ic,nL(ic)),*)Aatom(ic,nL(ic)) !remove space bef 
                  

                  if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU

                  if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B

                  if(du.eq.'_atom_site.auth_asym_id')i_ch=in !A,B, using later one

                  if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b

                  if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3

                  if(du.eq.'_atom_site.auth_seq_id')i_resi=in !1,2,3, identical to PDB res

                  if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?

                  

                  if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234
@@ -902,7 +918,8 @@ c read(Aatom(ic,nL(ic)),*)Aatom(ic,nL(ic)) !remove space bef 
                  endif

*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

                  

                  if(n_cut(ic).eq.0.and.nL(ic).gt.0)then !decide n_cut

                  if(n_cut(ic).eq.0.and.nL(ic).gt.0.and.

     &                 nLCA(ic).gt.0)then !decide n_cut for mmcif

                     read(ctmp(i_ch),*)ch_t

                     read(ctmp(i_ent),*)ent_t

                     if(ch_t.ne.Ach(ic,nL(ic)).or.ent_t.ne.
@@ -937,6 +954,9 @@ c^^^^^^^^^altLoc checked (mk.ne.1) will be skipped) ^^^^^^^^^ 
                     endif

c     read(ctmp(i_alt),*)Aalt(ic,nL(ic)) ! not used, because we alway use 1 atom for altLoc

c     if(Aalt(ic,nL(ic)).eq.'.')Aalt(ic,nL(ic))=' '

                     if(Aatom(ic,nL(ic)).eq.' CA')then

                        nLCA(ic)=nLCA(ic)+1 !for deciding n_cut of mmCIF

                     endif

                     

                     read(ctmp(i_res),*)Ares(ic,nL(ic))

                     read(ctmp(i_ch),*)Ach(ic,nL(ic)) !for check other chain
@@ -968,8 +988,6 @@ c if(Aalt(ic,nL(ic)).eq.'.')Aalt(ic,nL(ic))=' ' 
         close(10)

 1002 continue                  !ic=1,2 for file1 and file2

      

 8999 format(a6,I5,a1,A4,a1,A3,a1,A1,I4,A1,a3,3F8.3)



********** mark aligned residues --------------->

      do i=1,n8_al

         m12(1,i)=m1(i)
@@ -983,15 +1001,29 @@ c********************************************************************* 
c************* Step-2: output structures for Rasmol ******************

c*********************************************************************

      do 1019 ip=1,5

         if(ip.eq.1) outnameall=outname(1:len_trim(outname)) !'aaa', aligned CA

         if(ip.eq.2) outnameall=outname(1:len_trim(outname))//'_all' !'aaa_all', all CA

         if(ip.eq.3) outnameall=outname(1:len_trim(outname))//'_atm' !'aaa_atm', all aligned atoms

         if(ip.eq.4) outnameall=outname(1:len_trim(outname))

     &        //'_all_atm'      !'aaa_all_atm', all atoms

         if(ip.eq.5) outnameall=outname(1:len_trim(outname))

     &        //'_all_atm_lig'  !'aaa_all_atm_lig', all atom + ligand

         

         OPEN(unit=10,file=outnameall,status='unknown')

         if(ip.eq.1) xxx=outname(1:len_trim(outname)) !'aaa', aligned CA

         if(ip.eq.2) xxx=outname(1:len_trim(outname))//'_all' !'aaa_all', all CA

         if(ip.eq.3) xxx=outname(1:len_trim(outname))//'_atm' !'aaa_atm', all aligned atoms

         if(ip.eq.4) xxx=outname(1:len_trim(outname))//'_all_atm' !'aaa_all_atm', all atoms

         if(ip.eq.5) xxx=outname(1:len_trim(outname))//'_all_atm_lig' !'aaa_all_atm_lig', atom+ligand

         OPEN(unit=10,file=xxx,status='unknown')

         

*************for pymol preset ------>

         xxx_p=xxx(1:len_trim(xxx))//'.pml' !for pymol script

         xxx_pdb=xxx(1:len_trim(xxx))//'.pdb' !for pymol pdb

         OPEN(unit=11,file=xxx_p,status='unknown')

         OPEN(unit=12,file=xxx_pdb,status='unknown')

         

         write(11,'(A,A)')'load ',xxx_pdb(1:len_trim(xxx_pdb))

         write(11,'(A)')'hide all'

         write(11,'(A)')'bg_color white'

         write(11,'(A)')'color blue, chain A'

         write(11,'(A)')'color red, chain B'

         write(11,'(A)')'set transparency=0.2'

         write(11,'(A)')'set sphere_scale, 0.25'

         write(11,'(A)')'set stick_radius, 0.3'

*^^^^^^^^^ pymol preset complete ^^^^^^^^^^^^^^^^^^^

         

***   script:

         if(ip.eq.1.or.ip.eq.2)then !'aaa','aaa_all'

            write(10,'(A)')'load inline'
@@ -1010,6 +1042,9 @@ c********************************************************************* 
               endif

            enddo

            write(10,'(A)')'select all'

            

            write(11,'(A)')'show sticks, chain A'

            write(11,'(A)')'show sticks, chain B'

         elseif(ip.eq.3.or.ip.eq.4)then !'aaa_atm', 'aaa_all_atm'

            write(10,'(A)')'load inline'

            write(10,'(A)')'select *A'
@@ -1018,6 +1053,9 @@ c********************************************************************* 
            write(10,'(A)')'color red'

            write(10,'(A)')'select all'

            write(10,'(A)')'cartoon'

            

            write(11,'(A)')'show cartoon, chain A'

            write(11,'(A)')'show cartoon, chain B'

         else                   !'aaa_all_atm_lig

            write(10,'(A)')'load inline'

            write(10,'(A)')'select all'
@@ -1031,6 +1069,11 @@ c********************************************************************* 
            write(10,'(A)')'select solvent'

            write(10,'(A)')'spacefill 0.25'

            write(10,'(A)')'select all'

            

            write(11,'(A)')'show cartoon, chain A'

            write(11,'(A)')'show cartoon, chain B'

            write(11,'(A)')'show stick, HETATM'

            write(11,'(A)')'show spheres, solvent'

         endif

         write(10,'(A)')'exit'

         write(10,903)version
@@ -1045,11 +1088,11 @@ c------------output coordinates -----------> 
            else                !na=5001-10000

               na=5000          !for CONECT

            endif

c            write(*,*)'ip=',ip,'ic=',ic,nL(ic)

            

            if(n_cut(ic).eq.0)then

               n_cut(ic)=nL(ic) !there is no cut

            endif

            

            do i=1,nL(ic)

c---------- decide if we should output this line ------------>

               mk=0
@@ -1057,8 +1100,8 @@ c---------- decide if we should output this line ------------> 
                  if(i.le.n_cut(ic))then

                     if(Aatom(ic,i).eq.' CA')then

                        do j=1,n8_al

                           if(Aresi(ic,i).eq.mm(ic,m12(ic,j)))then

                              if(Ains(ic,i).eq.Cins(ic,m12(ic,j)))then

                           if(Aresi(ic,i).eq.mm(ic,m12(ic,j)))then !residue order

                              if(Ains(ic,i).eq.Cins(ic,m12(ic,j)))then !residue insertion

                                 mk=1

                                 goto 120 !each line output once

                              endif
@@ -1072,28 +1115,33 @@ c---------- decide if we should output this line ------------> 
                        mk=1

                     endif

                  endif

               elseif(ip.eq.3)then !'aaa_atm', all aligned atoms

               elseif(ip.eq.3)then !'aaa_atm', all aligned protein atoms

                  if(i.le.n_cut(ic))then

                     do j=1,n8_al

                        if(Aresi(ic,i).eq.mm(ic,m12(ic,j)))then

                           if(Ains(ic,i).eq.Cins(ic,m12(ic,j)))then

                        if(Ach(ic,i).eq.chid(ic,m12(ic,j)))then ! residue order (e.g. 100)

                           if(Aresi(ic,i).eq.mm(ic,m12(ic,j)))then ! residue order (e.g. 100)

                              if(Ains(ic,i).eq.Cins(ic,m12(ic,j)))then !resi insertion (eg, A)

                                 mk=1

                                 goto 120 !each line output once

                              endif

                           endif

                        endif

                     enddo

                  endif

               elseif(ip.eq.4)then !'aaa_all_atm', all atoms

               elseif(ip.eq.4)then !'aaa_all_atm', all protein atoms

                  if(Ach(ic,i).eq.chid(ic,1))then ! residue order (e.g. 100)

                     if(i.le.n_cut(ic))then

                        mk=1

                     endif

               elseif(ip.eq.5)then !'aaa_all_atm_lig', all atoms and ligands

                  endif

               elseif(ip.eq.5)then !'aaa_all_atm_lig', first chain and all ligands

                  if(i.le.n_cut(ic).or.Agroup(ic,i).eq.'HETATM')then

                     mk=1

                  endif

               endif

 120           continue



c^^^^^^^^^^

c               write(*,*)ip,ic,i,'mk=',mk,n_cut(ic)

               if(mk.eq.1)then  !mk=1----------->

                  if(ic.eq.1)then

                     ax=t(1)+u(1,1)*xx(ic,i)+u(1,2)*yy(ic,i)+
@@ -1113,29 +1161,40 @@ c--------------output (x,y,z) ---------------------------------> 
                  if(ip.eq.1.or.ip.eq.2)then !need CONECT

                     write(10,1236)Agroup(ic,i),na,'',

     &                    Aatom(ic,i),'',Ares(ic,i),'',Cch(ic),

     &                    Aresi(ic,i),Ains(ic,i),'',ax,ay,az

                     write(12,1236)Agroup(ic,i),na,'',

     &                    Aatom(ic,i),'',Ares(ic,i),'',Cch(ic),

     &                    Aresi(ic,i),Ains(ic,i),'',ax,ay,az

                  else          !all atoms

                     write(10,1236)Agroup(ic,i),Aatomi(ic,i),'',

     &                    Aatom(ic,i),'',Ares(ic,i),'',Cch(ic),

     &                    Aresi(ic,i),Ains(ic,i),'',ax,ay,az

                     write(12,1236)Agroup(ic,i),Aatomi(ic,i),'',

     &                    Aatom(ic,i),'',Ares(ic,i),'',Cch(ic),

     &                    Aresi(ic,i),Ains(ic,i),'',ax,ay,az

                  endif

               endif            !<----mk=1

            enddo               !do i=1,nL(ic)



            write(10,'(A)')'TER' !TER

            write(12,'(A)')'TER' !TER

            if(ip.eq.1.or.ip.eq.2)then

               if(ic.eq.1)then

                  do i=2,na

                     write(10,'(A6,I5,I5)')'CONECT',i-1,i !CONECT atom numbers

                     write(12,'(A6,I5,I5)')'CONECT',i-1,i !CONECT atom numbers

                  enddo

               else

                  do i=5002,na

                     write(10,'(A6,I5,I5)')'CONECT',i-1,i !CONECT atom numbers

                     write(12,'(A6,I5,I5)')'CONECT',i-1,i !CONECT atom numbers

                  enddo

               endif

            endif

         enddo                  !do ic=1,2

         close(10)

         close(11)

         close(12)

 1019 continue                  !ip=1,5



 902     format('select ',I4,':A,',I4,':B')

TMscore.cpp

0 → 100644
+5231 −0

File added.

Preview size limit exceeded, changes collapsed.

TMscore.f

+71 −11
Original line number Diff line number Diff line
@@ -40,6 +40,8 @@ 
*                 the output files.

*     2019/07/08: Enabled TM-score to support both PDB and mmCIF formats,

*                 and updated structure reading which makes program faster.

*     2019/08/18: Fixed multiple bugs associated with mmCIF formats.

*     2019/08/22: added output scripts for pymol, C++ version was included.

*************************************************************************

      

c        1         2         3         4         5         6         7 !
@@ -56,7 +58,7 @@ c 3456789012345678901234567890123456789012345678901234567890123456789012345678 
      common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score

      dimension k_ali(nmax),k_ali0(nmax)



      character*500 fnam,pdb(100),outname

      character*500 fnam,pdb(100),outname,xxx_p,xxx_pdb

      character*3 aa(-1:20),seqA(nmax),seqB(nmax),aanam,ss(2,nmax)

      character*500 s,du

      character*1 chA(nmax),chB(nmax),ch(nmax)
@@ -86,7 +88,7 @@ ccc mmCIF 
      character*20 Aatom(2,namax)

      character Agroup(2,namax)*6

      character Ares(2,namax)*3,Aalt(2,namax)

      character Ains(2,9000),Ach(2,namax),Aent(2,namax)

      character Ains(2,namax),Ach(2,namax),Aent(2,namax)

      character Cins(2,nmax)



      integer Aatomi(2,namax),Aresi(2,namax)
@@ -142,10 +144,13 @@ ccc 
         write(*,*)'5. TM-score with superposition output, e.g. ',

     &        '''TM.sup'':'

         write(*,*)'   >TMscore model native -o TM.sup'

         write(*,*)'   To view the superimposed CA-traces by rasmol:'

         write(*,*)'   To view superimposed CA-traces by rasmol ',

     &        'or pymol:'

         write(*,*)'      >rasmol -script TM.sup'

         write(*,*)'   To view superimposed atomic models by rasmol:'

         write(*,*)'      >pymol -d @TM.sup.pml'

         write(*,*)'   To view superimposed atomic models:'

         write(*,*)'      >rasmol -script TM.sup_atm'

         write(*,*)'      >pymol -d @TM.sup_atm.pml'

         write(*,*)

         goto 9999

      endif
@@ -284,8 +289,10 @@ c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^ 
                  if(du.eq.'_atom_site.label_alt_id')i_alt=in !'.',A,B

                  if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU

                  if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B

                  if(du.eq.'_atom_site.auth_asym_id')i_ch=in !A,B, using later one

                  if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b

                  if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3

                  if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3, 

                  if(du.eq.'_atom_site.auth_seq_id')i_resi=in !1,2,3, identical to PDB res

                  if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?

                  if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234

                  if(du.eq.'_atom_site.Cartn_y')i_y=in !y, 1.234
@@ -296,7 +303,7 @@ c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^ 
                  read(s,*)(ctmp(j),j=1,in) !no space before characters

                  if(ctmp(i_atom).eq.'CA')then !read 'CA' ---------->

                     if(i.gt.0)then

                        if(m_complex.eq.0)then

                        if(m_complex.eq.0)then !monomer

                           if(i_mn.gt.1)then !sometimes it may have no i_mn

                              read(ctmp(i_mn),*)mn_t

                              if(mn_t.ne.mn(i)) goto 1013 !only read first model
@@ -305,7 +312,7 @@ c^^^^^^^^^altLoc checked (mk.ne.1 will be skipped) ^^^^^^^^^ 
                           if(ch_t.ne.ch(i)) goto 1013 !only read first chain

                           read(ctmp(i_ent),*)ent_t

                           if(ent_t.ne.ent(i)) goto 1013 !only read first entity

                        else

                        else    !dimer

                           if(i_mn.gt.1)then !sometimes it may have no i_mn

                              read(ctmp(i_mn),*)mn_t

                              if(mn_t.ne.mn(i)) goto 1016 !only read first model
@@ -1026,13 +1033,13 @@ c^^^^^^^^^altLoc checked (mk.ne.1) will be skipped) ^^^^^^^^^ 
                  if(du.eq.'_atom_site.type_symbol')i_type=in !N,C,O

                  if(du.eq.'_atom_site.label_atom_id')i_atom=in !CA,O

                  if(du.eq.'_atom_site.label_alt_id')i_alt=in !'.',A,B

                  

                  if(du.eq.'_atom_site.label_comp_id')i_res=in !GLY,LEU

                  if(du.eq.'_atom_site.label_asym_id')i_ch=in !A,B

                  if(du.eq.'_atom_site.auth_asym_id')i_ch=in !A,B, using later one

                  if(du.eq.'_atom_site.label_entity_id')i_ent=in !1,2,a,b

                  if(du.eq.'_atom_site.label_seq_id')i_resi=in !1,2,3

                  if(du.eq.'_atom_site.auth_seq_id')i_resi=in !1,2,3, identical to PDB res

                  if(du.eq.'_atom_site.pdbx_PDB_ins_code')i_ins=in !A,?

                  

                  if(du.eq.'_atom_site.Cartn_x')i_x=in !x, 1.234

                  if(du.eq.'_atom_site.Cartn_y')i_y=in !y, 1.234

                  if(du.eq.'_atom_site.Cartn_z')i_z=in !z, 1.234
@@ -1122,6 +1129,24 @@ c************************************************************** 
c*************  Step-2: output 'aaa' (CA only)  ***************

c**************************************************************

      OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln



*************for pymol preset ------>

      xxx_p=outname(1:len_trim(outname))//'.pml' !for pymol script

      xxx_pdb=outname(1:len_trim(outname))//'.pdb' !for pymol script

      OPEN(unit=11,file=xxx_p,status='unknown')

      OPEN(unit=12,file=xxx_pdb,status='unknown')

      

      write(11,'(A,A)')'load ',xxx_pdb(1:len_trim(xxx_pdb))

      write(11,'(A)')'hide all'

      write(11,'(A)')'bg_color white'

      write(11,'(A)')'color blue, chain A'

      write(11,'(A)')'color red, chain B'

      write(11,'(A)')'set transparency=0.2'

      write(11,'(A)')'set stick_radius, 0.3'

      write(11,'(A)')'show sticks, chain A'

      write(11,'(A)')'show sticks, chain B'

*^^^^^^^^^ pymol preset complete ^^^^^^^^^^^^^^^^^^^

      

 900  format(A)

 901  format('select atomno= ',I5)

 902  format('select atomno= ',I5)
@@ -1146,21 +1171,31 @@ c************************************************************** 
 515  format('REMARK  TM-score=',f6.4,' (d0=',f5.2,')')

      do i=1,nseqA

         write(7,1236)i,seqA(i),chA(i),nresA(i),ins1(i),

     &        xt(i),yt(i),zt(i)

         write(12,1236)i,seqA(i),'A',nresA(i),ins1(i),

     &        xt(i),yt(i),zt(i)

      enddo

      write(7,1238)

      write(12,1238)

      do i=2,nseqA

         write(7,1239)i-1,i

         write(12,1239)i-1,i

      enddo

      do i=1,nseqB

         write(7,1236)50000+i,seqB(i),chB(i),nresB(i),ins2(i),

     &        xb(i),yb(i),zb(i)

         write(12,1236)50000+i,seqB(i),'B',nresB(i),ins2(i),

     &        xb(i),yb(i),zb(i)

      enddo

      write(7,1238)

      write(12,1238)

      do i=2,nseqB

         write(7,1239)50000+i-1,50000+i

         write(12,1239)50000+i-1,50000+i

      enddo

      close(7)

      close(11)

      close(12)

 1236 format('ATOM  ',i5,'  CA  ',A3,' ',A1,I4,A1,3X,3F8.3)

 1238 format('TER')

 1239 format('CONECT',I5,I5)
@@ -1169,9 +1204,26 @@ c************************************************************** 
c**************************************************************

c*************  Step-3: output 'aaa_atm' (all-atom) ***********

c**************************************************************

c     outname=trim(outname)//'_atm'

      outname=outname(1:len_trim(outname))//'_atm'

      OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln

      

*************for pymol preset ------>

      xxx_p=outname(1:len_trim(outname))//'.pml' !for pymol script

      xxx_pdb=outname(1:len_trim(outname))//'.pdb' !for pymol script

      OPEN(unit=11,file=xxx_p,status='unknown')

      OPEN(unit=12,file=xxx_pdb,status='unknown')

      

      write(11,'(A,A)')'load ',xxx_pdb(1:len_trim(xxx_pdb))

      write(11,'(A)')'hide all'

      write(11,'(A)')'bg_color white'

      write(11,'(A)')'color blue, chain A'

      write(11,'(A)')'color red, chain B'

      write(11,'(A)')'set transparency=0.2'

      write(11,'(A)')'set stick_radius, 0.3'

      write(11,'(A)')'show cartoon, chain A'

      write(11,'(A)')'show cartoon, chain B'

*^^^^^^^^^ pymol preset complete ^^^^^^^^^^^^^^^^^^^



      write(7,900)'load inline'

      write(7,900)'select atomno<50001'

      write(7,900)'color blue'
@@ -1189,16 +1241,24 @@ c outname=trim(outname)//'_atm' 
         az=t2(3)+u2(3,1)*xx(1,i)+u2(3,2)*yy(1,i)+u2(3,3)*zz(1,i)



         write(7,8888)i,Aatom(1,i),Ares(1,i),Ach(1,i),

     &        Aresi(1,i),Ains(1,i),ax,ay,az

         write(12,8888)i,Aatom(1,i),Ares(1,i),'A',

     &        Aresi(1,i),Ains(1,i),ax,ay,az

      enddo

      write(7,1238)             !TER

      write(12,1238)             !TER

***   chain2:

      do i=1,nL(2)

         write(7,8888)50000+i,Aatom(2,i),Ares(2,i),Ach(2,i),

     &        Aresi(2,i),Ains(2,i),xx(2,i),yy(2,i),zz(2,i)

         write(12,8888)50000+i,Aatom(2,i),Ares(2,i),'B',

     &        Aresi(2,i),Ains(2,i),xx(2,i),yy(2,i),zz(2,i)

      enddo

      write(7,1238)             !TER

      write(12,1238)             !TER

      close(7)

      close(11)

      close(12)

 8888 format('ATOM  ',I5,1x,A4,1x,A3,' ',A1,I4,A1,3x,3F8.3)

      

*^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^