rm uscdiff
[unres.git] / source / cluster / wham / src-M / readpdb.f
1       subroutine readpdb
2 C Read the PDB file and convert the peptide geometry into virtual-chain 
3 C geometry.
4       implicit real*8 (a-h,o-z)
5       include 'DIMENSIONS'
6       include 'COMMON.LOCAL'
7       include 'COMMON.VAR'
8       include 'COMMON.CHAIN'
9       include 'COMMON.INTERACT'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.GEO'
12       include 'COMMON.NAMES'
13       character*3 seq,atom,res
14       character*80 card
15       dimension sccor(3,20)
16       integer rescode
17 c      call permut(symetr)
18       ibeg=1
19       write(iout,*) 'pdbread'
20       do
21         read (ipdbin,'(a80)',end=10) card
22         if (card(:3).eq.'END') then
23           goto 10
24         else if (card(:3).eq.'TER') then
25 C End current chain
26           ires_old=ires+2
27           itype(ires_old-1)=ntyp1
28           itype(ires_old)=ntyp1
29           ibeg=2
30 c          write (iout,*) "Chain ended",ires,ishift,ires_old
31           call sccenter(ires,iii,sccor)
32         endif
33 C Fish out the ATOM cards.
34         if (index(card(1:4),'ATOM').gt.0) then  
35           read (card(14:16),'(a3)') atom
36           if (atom.eq.'CA' .or. atom.eq.'CH3') then
37 C Calculate the CM of the preceding residue.
38             if (ibeg.eq.0) then
39               call sccenter(ires,iii,sccor)
40             endif
41 C Start new residue.
42 c            write (iout,'(a80)') card
43             read (card(24:26),*) ires
44             read (card(18:20),'(a3)') res
45             if (ibeg.eq.1) then
46               ishift=ires-1
47               if (res.ne.'GLY' .and. res.ne. 'ACE') then
48                 ishift=ishift-1
49                 itype(1)=ntyp1
50               endif
51 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
52               ibeg=0          
53             else if (ibeg.eq.2) then
54 c Start a new chain
55               ishift=-ires_old+ires-1
56 c              write (iout,*) "New chain started",ires,ishift
57               ibeg=0
58             endif
59             ires=ires-ishift
60 c            write (2,*) "ires",ires," ishift",ishift
61             if (res.eq.'ACE') then
62               ity=10
63             else
64               itype(ires)=rescode(ires,res,0)
65             endif
66             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
67             write (iout,'(2i3,2x,a,3f8.3)') 
68      &       ires,itype(ires),res,(c(j,ires),j=1,3)
69             iii=1
70             do j=1,3
71               sccor(j,iii)=c(j,ires)
72             enddo
73           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
74      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
75             iii=iii+1
76             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
77           endif
78         endif
79       enddo
80    10 write (iout,'(a,i5)') ' Nres: ',ires
81 C Calculate dummy residue coordinates inside the "chain" of a multichain
82 C system
83       nres=ires
84       do i=2,nres-1
85 c        write (iout,*) i,itype(i)
86         if (itype(i).eq.ntyp1) then
87          if (itype(i+1).eq.ntyp1) then
88
89 c          write (iout,*) "dummy",i,itype(i)
90 C          do j=1,3
91 C            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
92 c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
93 C            dc(j,i)=c(j,i)
94 C          enddo
95            do j=1,3
96              dcj=(c(j,i-2)-c(j,i-3))/2.0
97              c(j,i)=c(j,i-1)+dcj
98              c(j,nres+i)=c(j,i)
99            enddo
100 C          endif   !unres_pdb
101          else     !itype(i+1).eq.ntyp1
102            do j=1,3
103             dcj=(c(j,i+3)-c(j,i+2))/2.0
104             c(j,i)=c(j,i+1)-dcj
105             c(j,nres+i)=c(j,i)
106            enddo
107 C          endif !unres_pdb
108          endif !itype(i+1).eq.ntyp1
109         endif  !itype.eq.ntyp1
110       enddo
111 C Calculate the CM of the last side chain.
112       call sccenter(ires,iii,sccor)
113       nsup=nres
114       nstart_sup=1
115       if (itype(nres).ne.10) then
116         nres=nres+1
117         itype(nres)=ntyp1
118         do j=1,3
119           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
120           c(j,nres)=c(j,nres-1)+dcj
121           c(j,2*nres)=c(j,nres)
122         enddo
123       endif
124       do i=2,nres-1
125         do j=1,3
126           c(j,i+nres)=dc(j,i)
127         enddo
128       enddo
129       do j=1,3
130         c(j,nres+1)=c(j,1)
131         c(j,2*nres)=c(j,nres)
132       enddo
133       if (itype(1).eq.ntyp1) then
134         nsup=nsup-1
135         nstart_sup=2
136         do j=1,3
137           dcj=(c(j,4)-c(j,3))/2.0
138           c(j,1)=c(j,2)-dcj
139           c(j,nres+1)=c(j,1)
140         enddo
141       endif
142 C Calculate internal coordinates.
143       do ires=1,nres
144         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') 
145      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
146      &    (c(j,nres+ires),j=1,3)
147       enddo
148       call int_from_cart1(.false.)
149       call int_from_cart(.true.,.false.)
150       call sc_loc_geom(.true.)
151       write (iout,*) "After int_from_cart"
152       call flush(iout)
153       do i=1,nres-1
154         do j=1,3
155           dc(j,i)=c(j,i+1)-c(j,i)
156           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
157         enddo
158       enddo
159       do i=2,nres-1
160         do j=1,3
161           dc(j,i+nres)=c(j,i+nres)-c(j,i)
162           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
163         enddo
164 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
165 c     &   vbld_inv(i+nres)
166       enddo
167       do i=1,nres
168         thetaref(i)=theta(i)
169         phiref(i)=phi(i)
170 c
171         phi_ref(i)=phi(i)
172         theta_ref(i)=theta(i)
173         alph_ref(i)=alph(i)
174         omeg_ref(i)=omeg(i)
175       enddo
176       
177 c      call chainbuild
178 C Copy the coordinates to reference coordinates
179 c      do i=1,2*nres
180 c        do j=1,3
181 c          cref(j,i)=c(j,i)
182 c        enddo
183 c      enddo
184
185       kkk=1
186       lll=0
187       cou=1
188       do i=1,nres
189       lll=lll+1
190 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
191       if (i.gt.1) then
192       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
193       chain_length=lll-1
194       kkk=kkk+1
195 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
196       lll=1
197       endif
198       endif
199         do j=1,3
200           cref(j,i,cou)=c(j,i)
201           cref(j,i+nres,cou)=c(j,i+nres)
202           if (i.le.nres) then
203           chain_rep(j,lll,kkk)=c(j,i)
204           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
205           endif
206          enddo
207       enddo
208       do j=1,3
209       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
210       chain_rep(j,chain_length+nres,symetr)
211      &=chain_rep(j,chain_length+nres,1)
212       enddo
213
214       if (symetr.gt.1) then
215        call permut(symetr)
216        nperm=1
217        do i=1,symetr
218        nperm=nperm*i
219        enddo
220 c       do i=1,nperm
221 c       write(iout,*) "tabperm", (tabperm(i,kkk),kkk=1,4)
222 c       enddo
223        do i=1,nperm
224         cou=0
225         do kkk=1,symetr
226          icha=tabperm(i,kkk)
227 c         write (iout,*) i,icha
228          do lll=1,chain_length
229           cou=cou+1
230            if (cou.le.nres) then
231            do j=1,3
232             kupa=mod(lll,chain_length)
233             iprzes=(kkk-1)*chain_length+lll
234             if (kupa.eq.0) kupa=chain_length
235 c            write (iout,*) "kupa", kupa
236             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
237             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
238           enddo
239           endif
240          enddo
241         enddo
242        enddo
243        endif
244
245 C-koniec robienia kopidm
246       nperm=0
247       do kkk=1,nperm
248       write (iout,*) "nowa struktura", nperm
249       do i=1,nres
250         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
251      &cref(2,i,kkk),
252      &cref(3,i,kkk),cref(1,nres+i,kkk),
253      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
254       enddo
255   100 format (//'              alpha-carbon coordinates       ',
256      &          '     centroid coordinates'/
257      1          '       ', 6X,'X',11X,'Y',11X,'Z',
258      &                          10X,'X',11X,'Y',11X,'Z')
259   110 format (a,'(',i3,')',6f12.5)
260        enddo
261
262
263       ishift_pdb=ishift
264       return
265       end
266 c---------------------------------------------------------------------------
267       subroutine int_from_cart(lside,lprn)
268       implicit real*8 (a-h,o-z)
269       include 'DIMENSIONS'
270       include 'COMMON.LOCAL'
271       include 'COMMON.VAR'
272       include 'COMMON.CHAIN'
273       include 'COMMON.INTERACT'
274       include 'COMMON.IOUNITS'
275       include 'COMMON.GEO'
276       include 'COMMON.NAMES'
277       character*3 seq,atom,res
278       character*80 card
279       dimension sccor(3,20)
280       integer rescode
281       logical lside,lprn
282       if (lprn) then 
283         write (iout,'(/a)') 
284      &  'Internal coordinates calculated from crystal structure.'
285         if (lside) then 
286           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
287      & '       Phi','    Dsc_id','       Dsc','     Alpha',
288      & '     Omega'
289         else 
290           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
291      & '       Phi'
292         endif
293       endif
294       call flush(iout)
295       do i=nnt+1,nct
296         iti=itype(i)
297 c        write (iout,*) i,dist(i,i-1)
298         if ((dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0)
299      &.and.(iti.ne.ntyp1).and.(itype(i-1).ne.ntyp1)) then
300           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
301           stop
302         endif
303         theta(i+1)=alpha(i-1,i,i+1)
304         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
305       enddo
306       if (lside) then
307         do i=2,nres-1
308           do j=1,3
309             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
310           enddo
311           iti=itype(i)
312           di=dist(i,nres+i)
313           if (iti.ne.10) then
314             alph(i)=alpha(nres+i,i,maxres2)
315             omeg(i)=beta(nres+i,i,maxres2,i+1)
316           endif
317           if (lprn)
318      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
319      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,rad2deg*alph(i),
320      &    rad2deg*omeg(i)
321         enddo
322       else if (lprn) then
323         do i=2,nres
324           iti=itype(i)
325           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
326      &    rad2deg*theta(i),rad2deg*phi(i)
327         enddo
328       endif
329       return
330       end
331
332 c-------------------------------------------------------------------------------
333       subroutine sc_loc_geom(lprn)
334       implicit real*8 (a-h,o-z)
335       include 'DIMENSIONS'
336       include 'COMMON.LOCAL'
337       include 'COMMON.VAR'
338       include 'COMMON.CHAIN'
339       include 'COMMON.INTERACT'
340       include 'COMMON.IOUNITS'
341       include 'COMMON.GEO'
342       include 'COMMON.NAMES'
343       include 'COMMON.CONTROL'
344       include 'COMMON.SETUP'
345       double precision x_prime(3),y_prime(3),z_prime(3)
346       logical lprn
347       do i=1,nres-1
348         do j=1,3
349           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
350         enddo
351       enddo
352       do i=2,nres-1
353         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
354           do j=1,3
355             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
356           enddo
357         else
358           do j=1,3
359             dc_norm(j,i+nres)=0.0d0
360           enddo
361         endif
362       enddo
363       do i=2,nres-1
364         costtab(i+1) =dcos(theta(i+1))
365         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
366         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
367         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
368         cosfac2=0.5d0/(1.0d0+costtab(i+1))
369         cosfac=dsqrt(cosfac2)
370         sinfac2=0.5d0/(1.0d0-costtab(i+1))
371         sinfac=dsqrt(sinfac2)
372         it=itype(i)
373         if (it.ne.10 .and. itype(i).ne.ntyp1) then
374 c
375 C  Compute the axes of tghe local cartesian coordinates system; store in
376 c   x_prime, y_prime and z_prime 
377 c
378         do j=1,3
379           x_prime(j) = 0.00
380           y_prime(j) = 0.00
381           z_prime(j) = 0.00
382         enddo
383         do j = 1,3
384           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
385           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
386         enddo
387         call vecpr(x_prime,y_prime,z_prime)
388 c
389 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
390 C to local coordinate system. Store in xx, yy, zz.
391 c
392         xx=0.0d0
393         yy=0.0d0
394         zz=0.0d0
395         do j = 1,3
396           xx = xx + x_prime(j)*dc_norm(j,i+nres)
397           yy = yy + y_prime(j)*dc_norm(j,i+nres)
398           zz = zz + z_prime(j)*dc_norm(j,i+nres)
399         enddo
400
401         xxref(i)=xx
402         yyref(i)=yy
403         zzref(i)=zz
404         else
405         xxref(i)=0.0d0
406         yyref(i)=0.0d0
407         zzref(i)=0.0d0
408         endif
409       enddo
410       if (lprn) then
411         do i=2,nres
412           iti=itype(i)
413           if(me.eq.king.or..not.out1file)
414      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
415      &      yyref(i),zzref(i)
416         enddo
417       endif
418       return
419       end
420 c---------------------------------------------------------------------------
421       subroutine sccenter(ires,nscat,sccor)
422       implicit real*8 (a-h,o-z)
423       include 'DIMENSIONS'
424       include 'COMMON.CHAIN'
425       dimension sccor(3,20)
426       do j=1,3
427         sccmj=0.0D0
428         do i=1,nscat
429           sccmj=sccmj+sccor(j,i) 
430         enddo
431         dc(j,ires)=sccmj/nscat
432       enddo
433       return
434       end
435 c----------------------------------------------------------------------
436       subroutine readpdb_template(k)
437 C Read the PDB file for read_constr_homology with read2sigma
438 C and convert the peptide geometry into virtual-chain geometry.
439       implicit real*8 (a-h,o-z)
440       include 'DIMENSIONS'
441       include 'COMMON.LOCAL'
442       include 'COMMON.VAR'
443       include 'COMMON.CHAIN'
444       include 'COMMON.INTERACT'
445       include 'COMMON.IOUNITS'
446       include 'COMMON.GEO'
447       include 'COMMON.NAMES'
448       include 'COMMON.CONTROL'
449       include 'COMMON.SETUP'
450       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
451       logical lprn /.false./,fail
452       double precision e1(3),e2(3),e3(3)
453       double precision dcj,efree_temp
454       character*3 seq,res
455       character*5 atom
456       character*80 card
457       double precision sccor(3,20)
458       integer rescode,iterter(maxres)
459       logical unres_pdb
460       unres_pdb=.false.
461       do i=1,maxres
462          iterter(i)=0
463       enddo
464       ibeg=1
465       ishift1=0
466       ishift=0
467 c      write (2,*) "UNRES_PDB",unres_pdb
468       ires=0
469       ires_old=0
470       iii=0
471       lsecondary=.false.
472       nhfrag=0
473       nbfrag=0
474       do
475         read (ipdbin,'(a80)',end=10) card
476         if (card(:3).eq.'END') then
477           goto 10
478         else if (card(:3).eq.'TER') then
479 C End current chain
480           ires_old=ires+2
481           itype(ires_old-1)=ntyp1 
482           iterter(ires_old-1)=1
483           itype(ires_old)=ntyp1
484           iterter(ires_old)=1
485           ibeg=2
486           write (iout,*) "Chain ended",ires,ishift,ires_old
487           if (unres_pdb) then
488             do j=1,3
489               dc(j,ires)=sccor(j,iii)
490             enddo
491           else 
492             call sccenter(ires,iii,sccor)
493           endif
494         endif
495 C Fish out the ATOM cards.
496         if (index(card(1:4),'ATOM').gt.0) then  
497           read (card(12:16),*) atom
498 c          write (iout,*) "! ",atom," !",ires
499 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
500           read (card(23:26),*) ires
501           read (card(18:20),'(a3)') res
502 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
503 c     &      " ires_old",ires_old
504 c          write (iout,*) "ishift",ishift," ishift1",ishift1
505 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
506           if (ires-ishift+ishift1.ne.ires_old) then
507 C Calculate the CM of the preceding residue.
508             if (ibeg.eq.0) then
509               if (unres_pdb) then
510                 do j=1,3
511                   dc(j,ires)=sccor(j,iii)
512                 enddo
513               else
514                 call sccenter(ires_old,iii,sccor)
515               endif
516               iii=0
517             endif
518 C Start new residue.
519             if (res.eq.'Cl-' .or. res.eq.'Na+') then
520               ires=ires_old
521               cycle
522             else if (ibeg.eq.1) then
523 c              write (iout,*) "BEG ires",ires
524               ishift=ires-1
525               if (res.ne.'GLY' .and. res.ne. 'ACE') then
526                 ishift=ishift-1
527                 itype(1)=ntyp1
528               endif
529               ires=ires-ishift+ishift1
530               ires_old=ires
531 c              write (iout,*) "ishift",ishift," ires",ires,
532 c     &         " ires_old",ires_old
533 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
534               ibeg=0          
535             else if (ibeg.eq.2) then
536 c Start a new chain
537               ishift=-ires_old+ires-1
538               ires=ires_old+1
539               write (iout,*) "New chain started",ires,ishift
540               ibeg=0          
541             else
542               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
543               ires=ires-ishift+ishift1
544               ires_old=ires
545             endif
546             if (res.eq.'ACE' .or. res.eq.'NHE') then
547               itype(ires)=10
548             else
549               itype(ires)=rescode(ires,res,0)
550             endif
551           else
552             ires=ires-ishift+ishift1
553           endif
554 c          write (iout,*) "ires_old",ires_old," ires",ires
555 c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
556 c            ishift1=ishift1+1
557 c          endif
558 c          write (2,*) "ires",ires," res ",res," ity",ity
559           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
560      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
561             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
562 c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
563 c#ifdef DEBUG
564 c            write (iout,'(2i3,2x,a,3f8.3)') 
565 c     &      ires,itype(ires),res,(c(j,ires),j=1,3)
566 c#endif
567             iii=iii+1
568             do j=1,3
569               sccor(j,iii)=c(j,ires)
570             enddo
571             if (ishift.ne.0) then
572               ires_ca=ires+ishift-ishift1
573             else
574               ires_ca=ires
575             endif
576 c            write (*,*) card(23:27),ires,itype(ires)
577           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
578      &             atom.ne.'N' .and. atom.ne.'C' .and.
579      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
580      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
581 c            write (iout,*) "sidechain ",atom
582             iii=iii+1
583             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
584           endif
585         endif
586       enddo
587    10 if(me.eq.king.or..not.out1file) 
588      & write (iout,'(a,i5)') ' Nres: ',ires
589 C Calculate dummy residue coordinates inside the "chain" of a multichain
590 C system
591       nres=ires
592       do i=2,nres-1
593 c        write (iout,*) i,itype(i),itype(i+1)
594         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
595          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
596 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
597 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
598 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
599            if (unres_pdb) then
600 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
601             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
602             if (fail) then
603               e2(1)=0.0d0
604               e2(2)=1.0d0
605               e2(3)=0.0d0
606             endif !fail
607             do j=1,3
608              c(j,i)=c(j,i-1)-1.9d0*e2(j)
609             enddo
610            else   !unres_pdb
611            do j=1,3
612              dcj=(c(j,i-2)-c(j,i-3))/2.0
613             if (dcj.eq.0) dcj=1.23591524223
614              c(j,i)=c(j,i-1)+dcj
615              c(j,nres+i)=c(j,i)
616            enddo     
617           endif   !unres_pdb
618          else     !itype(i+1).eq.ntyp1
619           if (unres_pdb) then
620 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
621             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
622             if (fail) then
623               e2(1)=0.0d0
624               e2(2)=1.0d0
625               e2(3)=0.0d0
626             endif
627             do j=1,3
628               c(j,i)=c(j,i+1)-1.9d0*e2(j)
629             enddo
630           else !unres_pdb
631            do j=1,3
632             dcj=(c(j,i+3)-c(j,i+2))/2.0
633             if (dcj.eq.0) dcj=1.23591524223
634             c(j,i)=c(j,i+1)-dcj
635             c(j,nres+i)=c(j,i)
636            enddo
637           endif !unres_pdb
638          endif !itype(i+1).eq.ntyp1
639         endif  !itype.eq.ntyp1
640       enddo
641 C Calculate the CM of the last side chain.
642       if (unres_pdb) then
643         do j=1,3
644           dc(j,ires)=sccor(j,iii)
645         enddo
646       else
647         call sccenter(ires,iii,sccor)
648       endif
649       nsup=nres
650       nstart_sup=1
651       if (itype(nres).ne.10) then
652         nres=nres+1
653         itype(nres)=ntyp1
654         if (unres_pdb) then
655 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
656           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
657           if (fail) then
658             e2(1)=0.0d0
659             e2(2)=1.0d0
660             e2(3)=0.0d0
661           endif
662           do j=1,3
663             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
664           enddo
665         else
666         do j=1,3
667           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
668             if (dcj.eq.0) dcj=1.23591524223
669           c(j,nres)=c(j,nres-1)+dcj
670           c(j,2*nres)=c(j,nres)
671         enddo
672       endif
673       endif
674       do i=2,nres-1
675         do j=1,3
676           c(j,i+nres)=dc(j,i)
677         enddo
678       enddo
679       do j=1,3
680         c(j,nres+1)=c(j,1)
681         c(j,2*nres)=c(j,nres)
682       enddo
683       if (itype(1).eq.ntyp1) then
684         nsup=nsup-1
685         nstart_sup=2
686         if (unres_pdb) then
687 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
688           call refsys(2,3,4,e1,e2,e3,fail)
689           if (fail) then
690             e2(1)=0.0d0
691             e2(2)=1.0d0
692             e2(3)=0.0d0
693           endif
694           do j=1,3
695             c(j,1)=c(j,2)-1.9d0*e2(j)
696           enddo
697         else
698         do j=1,3
699           dcj=(c(j,4)-c(j,3))/2.0
700           c(j,1)=c(j,2)-dcj
701           c(j,nres+1)=c(j,1)
702         enddo
703         endif
704       endif
705 C Copy the coordinates to reference coordinates
706 c      do i=1,2*nres
707 c        do j=1,3
708 c          cref(j,i)=c(j,i)
709 c        enddo
710 c      enddo
711 C Calculate internal coordinates.
712       if (lprn) then
713       write (iout,'(/a)') 
714      &  "Cartesian coordinates of the reference structure"
715       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
716      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
717       do ires=1,nres
718         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') 
719      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
720      &    (c(j,ires+nres),j=1,3)
721       enddo
722       endif
723 C Calculate internal coordinates.
724       if(me.eq.king.or..not.out1file)then
725        write (iout,'(a)') 
726      &   "Backbone and SC coordinates as read from the PDB"
727        do ires=1,nres
728         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') 
729      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
730      &    (c(j,nres+ires),j=1,3)
731        enddo
732       endif
733       call int_from_cart1(.false.)
734       call int_from_cart(.true.,.false.)
735       call sc_loc_geom(.false.)
736       do i=1,nres
737         thetaref(i)=theta(i)
738         phiref(i)=phi(i)
739       enddo
740       do i=1,nres-1
741         do j=1,3
742           dc(j,i)=c(j,i+1)-c(j,i)
743           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
744         enddo
745       enddo
746       do i=2,nres-1
747         do j=1,3
748           dc(j,i+nres)=c(j,i+nres)-c(j,i)
749           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
750         enddo
751 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
752 c     &   vbld_inv(i+nres)
753       enddo
754 c      call chainbuild
755 C Copy the coordinates to reference coordinates
756 C Splits to single chain if occurs
757       kkk=1
758       lll=0
759       cou=1
760       do i=1,nres
761       lll=lll+1
762 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
763       if (i.gt.1) then
764       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
765       chain_length=lll-1
766       kkk=kkk+1
767 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
768       lll=1
769       endif
770       endif
771         do j=1,3
772           cref(j,i,cou)=c(j,i)
773           cref(j,i+nres,cou)=c(j,i+nres)
774           if (i.le.nres) then
775           chain_rep(j,lll,kkk)=c(j,i)
776           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
777           endif
778         enddo
779       enddo
780 c      do i=1,2*nres
781 c        do j=1,3
782 c          chomo(j,i,k)=c(j,i)
783 c        enddo
784 c      enddo
785
786       write (iout,*) chain_length
787       if (chain_length.eq.0) chain_length=nres
788       do j=1,3
789       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
790       chain_rep(j,chain_length+nres,symetr)
791      &=chain_rep(j,chain_length+nres,1)
792       enddo
793 c diagnostic
794 c       write (iout,*) "spraw lancuchy",chain_length,symetr
795 c       do i=1,4
796 c         do kkk=1,chain_length
797 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
798 c         enddo
799 c        enddo
800 c enddiagnostic       
801
802       return
803       end