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