5 implicit real*8 (a-h,o-z)
7 include 'DIMENSIONS.ZSCOPT'
8 include 'DIMENSIONS.FREE'
9 include 'COMMON.IOUNITS'
12 include 'COMMON.INTERACT'
13 include 'COMMON.LOCAL'
14 include 'COMMON.NAMES'
15 include 'COMMON.CHAIN'
16 include 'COMMON.FFIELD'
17 include 'COMMON.SBRIDGE'
18 include 'COMMON.TORCNSTR'
19 include 'COMMON.CONTROL'
20 character*4 sequence(maxres)
22 double precision x(maxvar)
23 character*320 controlcard,ucase
24 dimension itype_pdb(maxres)
26 call card_concat(controlcard,.true.)
27 call reada(controlcard,'SCAL14',scal14,0.4d0)
28 call reada(controlcard,'SCALSCP',scalscp,1.0d0)
29 call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
30 call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
31 r0_corr=cutoff_corr-delt_corr
32 call readi(controlcard,"NRES",nres,0)
33 iscode=index(controlcard,"ONE_LETTER")
35 write (iout,*) "Error: no residues in molecule"
38 if (nres.gt.maxres) then
39 write (iout,*) "Error: too many residues",nres,maxres
41 write(iout,*) 'nres=',nres
42 C Read sequence of the protein
44 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
46 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
48 C Convert sequence to numeric code
50 itype(i)=rescode(i,sequence(i),iscode)
52 write (iout,*) "Numeric code:"
53 write (iout,'(20i4)') (itype(i),i=1,nres)
56 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
58 if (itype(i).eq.ntyp1) then
62 else if (iabs(itype(i+1)).ne.20) then
64 else if (iabs(itype(i)).ne.20) then
73 write (iout,*) i,itype(i),itel(i)
77 if (with_dihed_constr) then
79 read (inp,*) ndih_constr
80 if (ndih_constr.gt.0) then
82 C write (iout,*) 'FTORS',ftors
83 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
86 & 'There are',ndih_constr,' constraints on phi angles.'
88 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
92 phi0(i)=deg2rad*phi0(i)
93 drange(i)=deg2rad*drange(i)
98 if (with_theta_constr) then
99 C with_theta_constr is keyword allowing for occurance of theta constrains
100 read (inp,*) ntheta_constr
101 C ntheta_constr is the number of theta constrains
102 if (ntheta_constr.gt.0) then
104 read (inp,*) (itheta_constr(i),theta_constr0(i),
105 & theta_drange(i),for_thet_constr(i),
107 C the above code reads from 1 to ntheta_constr
108 C itheta_constr(i) residue i for which is theta_constr
109 C theta_constr0 the global minimum value
110 C theta_drange is range for which there is no energy penalty
111 C for_thet_constr is the force constant for quartic energy penalty
113 C if(me.eq.king.or..not.out1file)then
115 & 'There are',ntheta_constr,' constraints on phi angles.'
117 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
123 theta_constr0(i)=deg2rad*theta_constr0(i)
124 theta_drange(i)=deg2rad*theta_drange(i)
126 C if(me.eq.king.or..not.out1file)
127 C & write (iout,*) 'FTORS',ftors
128 C do i=1,ntheta_constr
129 C ii = itheta_constr(i)
130 C thetabound(1,ii) = phi0(i)-drange(i)
131 C thetabound(2,ii) = phi0(i)+drange(i)
133 endif ! ntheta_constr.gt.0
134 endif! with_theta_constr
137 if (itype(1).eq.ntyp1) nnt=2
138 if (itype(nres).eq.ntyp1) nct=nct-1
139 write(iout,*) 'NNT=',NNT,' NCT=',NCT
140 write (iout,*) "calling read_saxs_consrtr",nsaxs
141 if (nsaxs.gt.0) call read_saxs_constr
143 if (constr_homology.gt.0) then
144 c write (iout,*) "About to call read_constr_homology"
146 call read_constr_homology
147 write (iout,*) "Exit read_constr_homology"
149 cref if (indpdb.gt.0 .or. pdbref) then
152 cref c(j,i)=crefjlee(j,i)
153 cref cref(j,i)=crefjlee(j,i)
165 write (iout,'(/a,i3,a)') 'The chain contains',ns,
166 & ' disulfide-bridging cysteines.'
167 write (iout,'(20i4)') (iss(i),i=1,ns)
169 write(iout,*)"Running with dynamic disulfide-bond formation"
171 write (iout,'(/a/)') 'Pre-formed links are:'
177 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
178 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',
179 & dhpb(i),ebr,forcon(i)
184 if (ns.gt.0.and.dyn_ss) then
188 forcon(i-nss)=forcon(i)
195 dyn_ss_mask(iss(i))=.true.
200 c-----------------------------------------------------------------------------
201 logical function seq_comp(itypea,itypeb,length)
203 integer length,itypea(length),itypeb(length)
206 if (itypea(i).ne.itypeb(i)) then
214 c-----------------------------------------------------------------------------
215 subroutine read_bridge
216 C Read information about disulfide bridges.
217 implicit real*8 (a-h,o-z)
219 include 'DIMENSIONS.ZSCOPT'
220 include 'COMMON.IOUNITS'
223 include 'COMMON.INTERACT'
224 include 'COMMON.NAMES'
225 include 'COMMON.CHAIN'
226 include 'COMMON.FFIELD'
227 include 'COMMON.SBRIDGE'
228 C Read bridging residues.
229 read (inp,*) ns,(iss(i),i=1,ns)
231 write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
232 C Check whether the specified bridging residues are cystines.
234 if (itype(iss(i)).ne.1) then
235 write (iout,'(2a,i3,a)')
236 & 'Do you REALLY think that the residue ',
237 & restyp(itype(iss(i))),i,
238 & ' can form a disulfide bridge?!!!'
239 write (*,'(2a,i3,a)')
240 & 'Do you REALLY think that the residue ',
241 & restyp(itype(iss(i))),i,
242 & ' can form a disulfide bridge?!!!'
246 C Read preformed bridges.
248 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
249 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
252 C Check if the residues involved in bridges are in the specified list of
256 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
257 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
258 write (iout,'(a,i3,a)') 'Disulfide pair',i,
259 & ' contains residues present in other pairs.'
260 write (*,'(a,i3,a)') 'Disulfide pair',i,
261 & ' contains residues present in other pairs.'
266 if (ihpb(i).eq.iss(j)) goto 10
268 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
271 if (jhpb(i).eq.iss(j)) goto 20
273 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
286 c------------------------------------------------------------------------------
287 subroutine read_angles(kanal,iscor,energ,iprot,*)
288 implicit real*8 (a-h,o-z)
290 include 'DIMENSIONS.ZSCOPT'
291 include 'COMMON.INTERACT'
292 include 'COMMON.SBRIDGE'
295 include 'COMMON.CHAIN'
296 include 'COMMON.IOUNITS'
298 read(kanal,'(a80)',end=10,err=10) lineh
299 read(lineh(:5),*,err=8) ic
300 read(lineh(6:),*,err=8) energ
303 print *,'error, assuming e=1d10',lineh
307 read(lineh(18:),*,end=10,err=10) nss
309 read (lineh(20:),*,end=10,err=10)
310 & (IHPB(I),JHPB(I),I=1,NSS),iscor
312 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
313 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),
316 c print *,"energy",energ," iscor",iscor
317 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
318 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
319 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
320 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
322 theta(i)=deg2rad*theta(i)
323 phi(i)=deg2rad*phi(i)
324 alph(i)=deg2rad*alph(i)
325 omeg(i)=deg2rad*omeg(i)
330 c-------------------------------------------------------------------------------
331 subroutine read_saxs_constr
332 implicit real*8 (a-h,o-z)
334 include 'DIMENSIONS.ZSCOPT'
335 include 'DIMENSIONS.FREE'
339 include 'COMMON.CONTROL'
340 include 'COMMON.CHAIN'
341 include 'COMMON.IOUNITS'
342 include 'COMMON.SBRIDGE'
343 double precision cm(3)
345 write (iout,*) "Calling read_saxs nsaxs",nsaxs
347 if (saxs_mode.eq.0) then
348 c SAXS distance distribution
350 read(inp,*) distsaxs(i),Psaxs(i)
354 Cnorm = Cnorm + Psaxs(i)
356 write (iout,*) "Cnorm",Cnorm
358 Psaxs(i)=Psaxs(i)/Cnorm
360 write (iout,*) "Normalized distance distribution from SAXS"
362 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
366 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
368 write (iout,*) "Wsaxs0",Wsaxs0
372 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
379 cm(j)=cm(j)+Csaxs(j,i)
387 Csaxs(j,i)=Csaxs(j,i)-cm(j)
390 write (iout,*) "SAXS sphere coordinates"
392 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
397 c====-------------------------------------------------------------------
398 subroutine read_constr_homology
401 include 'DIMENSIONS.ZSCOPT'
402 include 'DIMENSIONS.FREE'
406 include 'COMMON.SETUP'
407 include 'COMMON.CONTROL'
408 include 'COMMON.CHAIN'
409 include 'COMMON.IOUNITS'
411 include 'COMMON.INTERACT'
412 include 'include_unres/COMMON.NAMES'
413 include 'COMMON.HOMRESTR'
418 c include 'include_unres/COMMON.VAR'
421 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
423 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
424 c & sigma_odl_temp(maxres,maxres,max_template)
426 character*24 model_ki_dist, model_ki_angle
427 character*500 controlcard
428 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
429 integer idomain(max_template,maxres)
430 logical lprn /.true./
433 logical unres_pdb,liiflag
435 c FP - Nov. 2014 Temporary specifications for new vars
437 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
439 double precision, dimension (max_template,maxres) :: rescore
440 double precision, dimension (max_template,maxres) :: rescore2
441 double precision, dimension (max_template,maxres) :: rescore3
442 character*24 tpl_k_rescore
443 c -----------------------------------------------------------------
444 c Reading multiple PDB ref structures and calculation of retraints
445 c not using pre-computed ones stored in files model_ki_{dist,angle}
447 c -----------------------------------------------------------------
450 c Alternative: reading from input
451 call card_concat(controlcard,.true.)
452 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
453 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
454 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
455 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
456 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
457 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
458 dist1cut=(index(controlcard,'DIST1CUT').gt.0)
459 call readi(controlcard,"HOMOL_SET",homol_nset,1)
460 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
461 call readi(controlcard,"IHSET",ihset,1)
462 if (homol_nset.gt.1)then
463 call card_concat(controlcard,.true.)
464 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
465 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
466 write(iout,*) "iset homology_weight "
468 c write(iout,*) i,waga_homology(i)
471 iset=mod(kolor,homol_nset)+1
476 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
478 cd write (iout,*) "nnt",nnt," nct",nct
488 c Reading HM global scores (prob not required)
491 do k=1,constr_homology
495 c open (4,file="HMscore")
496 c do k=1,constr_homology
497 c read (4,*,end=521) hmscore_tmp
498 c hmscore(k)=hmscore_tmp ! Another transformation can be used
499 c write(*,*) "Model", k, ":", hmscore(k)
510 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
512 if (read_homol_frag) then
513 call read_klapaucjusz
516 do k=1,constr_homology
518 read(inp,'(a)') pdbfile
519 c Next stament causes error upon compilation (?)
520 c if(me.eq.king.or. .not. out1file)
521 c write (iout,'(2a)') 'PDB data will be read from file ',
522 c & pdbfile(:ilen(pdbfile))
523 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
524 & pdbfile(:ilen(pdbfile))
525 open(ipdbin,file=pdbfile,status='old',err=33)
527 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
528 & pdbfile(:ilen(pdbfile))
531 c print *,'Begin reading pdb data'
533 c Files containing res sim or local scores (former containing sigmas)
536 write(kic2,'(bz,i2.2)') k
538 tpl_k_rescore="template"//kic2//".sco"
541 call readpdb_template(k)
544 cref crefjlee(j,i)=c(j,i)
549 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
550 & (crefjlee(j,i+nres),j=1,3)
553 write (iout,*) "read_constr_homology: after reading pdb file"
557 c Distance restraints
560 C Copy the coordinates from reference coordinates (?)
564 c write (iout,*) "c(",j,i,") =",c(j,i)
568 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
570 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
571 open (ientin,file=tpl_k_rescore,status='old')
572 if (nnt.gt.1) rescore(k,1)=0.0d0
573 do irec=nnt,nct ! loop for reading res sim
575 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
576 & rescore3_tmp,idomain_tmp
578 idomain(k,i_tmp)=idomain_tmp
579 rescore(k,i_tmp)=rescore_tmp
580 rescore2(k,i_tmp)=rescore2_tmp
581 rescore3(k,i_tmp)=rescore3_tmp
582 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
583 & i_tmp,rescore2_tmp,rescore_tmp,
584 & rescore3_tmp,idomain_tmp
587 read (ientin,*,end=1401) rescore_tmp
589 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
590 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
591 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
596 if (waga_dist.ne.0.0d0) then
604 distal=dsqrt(x12*x12+y12*y12+z12*z12)
605 c write (iout,*) k,i,j,distal,dist2_cut
606 if (dist1cut .and. k.gt.1) then
608 if (l_homo(1,ii)) then
614 sigma_odl(k,ii)=sigma_odl(1,ii)
619 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
620 & .and. distal.le.dist2_cut ) then
626 c write (iout,*) "k",k
627 c write (iout,*) "i",i," j",j," constr_homology",
635 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
637 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
638 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
639 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
641 if (odl(k,ii).le.dist_cut) then
642 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
645 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
646 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
648 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
649 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
653 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
664 c Theta, dihedral and SC retraints
666 if (waga_angle.gt.0.0d0) then
667 c open (ientin,file=tpl_k_sigma_dih,status='old')
668 c do irec=1,maxres-3 ! loop for reading sigma_dih
669 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
670 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
671 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
672 c & sigma_dih(k,i+nnt-1)
677 if (idomain(k,i).eq.0) then
681 dih(k,i)=phiref(i) ! right?
682 c read (ientin,*) sigma_dih(k,i) ! original variant
683 c write (iout,*) "dih(",k,i,") =",dih(k,i)
684 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
685 c & "rescore(",k,i-1,") =",rescore(k,i-1),
686 c & "rescore(",k,i-2,") =",rescore(k,i-2),
687 c & "rescore(",k,i-3,") =",rescore(k,i-3)
689 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
690 & rescore(k,i-2)+rescore(k,i-3))/4.0
691 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
692 c write (iout,*) "Raw sigmas for dihedral angle restraints"
693 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
694 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
695 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
696 c Instead of res sim other local measure of b/b str reliability possible
697 if (sigma_dih(k,i).ne.0)
698 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
699 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
704 if (waga_theta.gt.0.0d0) then
705 c open (ientin,file=tpl_k_sigma_theta,status='old')
706 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
707 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
708 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
709 c & sigma_theta(k,i+nnt-1)
714 do i = nnt+2,nct ! right? without parallel.
715 c do i = i=1,nres ! alternative for bounds acc to readpdb?
716 c do i=ithet_start,ithet_end ! with FG parallel.
717 if (idomain(k,i).eq.0) then
721 thetatpl(k,i)=thetaref(i)
722 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
723 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
724 c & "rescore(",k,i-1,") =",rescore(k,i-1),
725 c & "rescore(",k,i-2,") =",rescore(k,i-2)
726 c read (ientin,*) sigma_theta(k,i) ! 1st variant
727 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
728 & rescore(k,i-2))/3.0
729 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
730 if (sigma_theta(k,i).ne.0)
731 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
733 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
734 c rescore(k,i-2) ! right expression ?
735 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
739 if (waga_d.gt.0.0d0) then
740 c open (ientin,file=tpl_k_sigma_d,status='old')
741 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
742 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
743 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
744 c & sigma_d(k,i+nnt-1)
748 do i = nnt,nct ! right? without parallel.
749 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
750 c do i=loc_start,loc_end ! with FG parallel.
751 if (itype(i).eq.10) cycle
752 if (idomain(k,i).eq.0 ) then
759 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
760 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
761 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
762 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
763 sigma_d(k,i)=rescore3(k,i) ! right expression ?
764 if (sigma_d(k,i).ne.0)
765 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
767 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
768 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
769 c read (ientin,*) sigma_d(k,i) ! 1st variant
774 c remove distance restraints not used in any model from the list
775 c shift data in all arrays
777 if (waga_dist.ne.0.0d0) then
783 if (ii_in_use(ii).eq.0.and.liiflag) then
787 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
788 & .not.liiflag.and.ii.eq.lim_odl) then
789 if (ii.eq.lim_odl) then
795 do ki=iistart,lim_odl-iishift
796 ires_homo(ki)=ires_homo(ki+iishift)
797 jres_homo(ki)=jres_homo(ki+iishift)
798 ii_in_use(ki)=ii_in_use(ki+iishift)
799 do k=1,constr_homology
800 odl(k,ki)=odl(k,ki+iishift)
801 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
802 l_homo(k,ki)=l_homo(k,ki+iishift)
806 lim_odl=lim_odl-iishift
812 endif ! .not. klapaucjusz
814 if (constr_homology.gt.0) call homology_partition
815 if (constr_homology.gt.0) call init_int_table
816 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
817 cd & "lim_xx=",lim_xx
818 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
819 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
823 if (.not.lprn) return
824 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
825 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
826 write (iout,*) "Distance restraints from templates"
828 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
829 & ii,ires_homo(ii),jres_homo(ii),
830 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
831 & ki=1,constr_homology)
833 write (iout,*) "Dihedral angle restraints from templates"
835 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
836 & (rad2deg*dih(ki,i),
837 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
839 write (iout,*) "Virtual-bond angle restraints from templates"
841 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
842 & (rad2deg*thetatpl(ki,i),
843 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
845 write (iout,*) "SC restraints from templates"
847 write(iout,'(i5,100(4f8.2,4x))') i,
848 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
849 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
852 c -----------------------------------------------------------------
855 c----------------------------------------------------------------------
856 subroutine read_klapaucjusz
859 include 'DIMENSIONS.ZSCOPT'
860 include 'DIMENSIONS.FREE'
864 include 'COMMON.SETUP'
865 include 'COMMON.CONTROL'
866 include 'COMMON.CHAIN'
867 include 'COMMON.IOUNITS'
869 include 'COMMON.INTERACT'
870 include 'COMMON.NAMES'
871 include 'COMMON.HOMRESTR'
872 character*256 fragfile
873 integer ninclust(maxclust),inclust(max_template,maxclust),
874 & nresclust(maxclust),iresclust(maxres,maxclust)
877 character*24 model_ki_dist, model_ki_angle
878 character*500 controlcard
879 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
880 integer idomain(max_template,maxres)
881 logical lprn /.true./
884 logical unres_pdb,liiflag
887 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
888 double precision, dimension (max_template,maxres) :: rescore
889 double precision, dimension (max_template,maxres) :: rescore2
890 character*24 tpl_k_rescore
897 double precision chomo(3,maxres2+2,max_template)
898 call getenv("FRAGFILE",fragfile)
899 open(ientin,file=fragfile,status="old",err=10)
900 read(ientin,*) constr_homology,nclust
906 do k=1,constr_homology
907 read(ientin,'(a)') pdbfile
908 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
909 & pdbfile(:ilen(pdbfile))
910 open(ipdbin,file=pdbfile,status='old',err=33)
912 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
913 & pdbfile(:ilen(pdbfile))
917 call readpdb_template(k)
930 read(ientin,*) ninclust(i),nresclust(i)
931 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
932 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
938 do ll = 1,ninclust(l)
946 idomain(k,iresclust(i,l)+1) = 1
948 idomain(k,iresclust(i,l)) = 1
952 c Distance restraints
955 C Copy the coordinates from reference coordinates (?)
959 c write (iout,*) "c(",j,i,") =",c(j,i)
962 call int_from_cart1(.false.)
963 call int_from_cart(.true.,.false.)
964 call sc_loc_geom(.false.)
969 if (waga_dist.ne.0.0d0) then
977 distal=dsqrt(x12*x12+y12*y12+z12*z12)
978 c write (iout,*) k,i,j,distal,dist2_cut
980 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
981 & .and. distal.le.dist2_cut ) then
987 c write (iout,*) "k",k
988 c write (iout,*) "i",i," j",j," constr_homology",
996 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
998 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
999 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
1000 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1002 if (odl(k,ii).le.dist_cut) then
1003 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1006 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1007 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1009 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
1010 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1014 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1017 c l_homo(k,ii)=.false.
1024 c Theta, dihedral and SC retraints
1026 if (waga_angle.gt.0.0d0) then
1028 if (idomain(k,i).eq.0) then
1029 c sigma_dih(k,i)=0.0
1033 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1034 & rescore(k,i-2)+rescore(k,i-3))/4.0
1035 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1036 c & " sigma_dihed",sigma_dih(k,i)
1037 if (sigma_dih(k,i).ne.0)
1038 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1043 if (waga_theta.gt.0.0d0) then
1045 if (idomain(k,i).eq.0) then
1046 c sigma_theta(k,i)=0.0
1049 thetatpl(k,i)=thetaref(i)
1050 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1051 & rescore(k,i-2))/3.0
1052 if (sigma_theta(k,i).ne.0)
1053 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1057 if (waga_d.gt.0.0d0) then
1059 if (itype(i).eq.10) cycle
1060 if (idomain(k,i).eq.0 ) then
1067 sigma_d(k,i)=rescore(k,i)
1068 if (sigma_d(k,i).ne.0)
1069 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1070 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1076 c remove distance restraints not used in any model from the list
1077 c shift data in all arrays
1079 if (waga_dist.ne.0.0d0) then
1085 if (ii_in_use(ii).eq.0.and.liiflag) then
1089 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1090 & .not.liiflag.and.ii.eq.lim_odl) then
1091 if (ii.eq.lim_odl) then
1092 iishift=ii-iistart+1
1097 do ki=iistart,lim_odl-iishift
1098 ires_homo(ki)=ires_homo(ki+iishift)
1099 jres_homo(ki)=jres_homo(ki+iishift)
1100 ii_in_use(ki)=ii_in_use(ki+iishift)
1101 do k=1,constr_homology
1102 odl(k,ki)=odl(k,ki+iishift)
1103 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1104 l_homo(k,ki)=l_homo(k,ki+iishift)
1108 lim_odl=lim_odl-iishift
1115 10 stop "Error infragment file"
1117 c----------------------------------------------------------------------