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 call readi(controlcard,"HOMOL_SET",homol_nset,1)
459 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
460 call readi(controlcard,"IHSET",ihset,1)
461 if (homol_nset.gt.1)then
462 call card_concat(controlcard,.true.)
463 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
464 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
465 write(iout,*) "iset homology_weight "
467 c write(iout,*) i,waga_homology(i)
470 iset=mod(kolor,homol_nset)+1
475 c write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
477 cd write (iout,*) "nnt",nnt," nct",nct
487 c Reading HM global scores (prob not required)
490 do k=1,constr_homology
494 c open (4,file="HMscore")
495 c do k=1,constr_homology
496 c read (4,*,end=521) hmscore_tmp
497 c hmscore(k)=hmscore_tmp ! Another transformation can be used
498 c write(*,*) "Model", k, ":", hmscore(k)
509 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
511 if (read_homol_frag) then
512 call read_klapaucjusz
515 do k=1,constr_homology
517 read(inp,'(a)') pdbfile
518 c Next stament causes error upon compilation (?)
519 c if(me.eq.king.or. .not. out1file)
520 c write (iout,'(2a)') 'PDB data will be read from file ',
521 c & pdbfile(:ilen(pdbfile))
522 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
523 & pdbfile(:ilen(pdbfile))
524 open(ipdbin,file=pdbfile,status='old',err=33)
526 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
527 & pdbfile(:ilen(pdbfile))
530 c print *,'Begin reading pdb data'
532 c Files containing res sim or local scores (former containing sigmas)
535 write(kic2,'(bz,i2.2)') k
537 tpl_k_rescore="template"//kic2//".sco"
540 call readpdb_template(k)
543 cref crefjlee(j,i)=c(j,i)
548 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
549 & (crefjlee(j,i+nres),j=1,3)
552 write (iout,*) "read_constr_homology: after reading pdb file"
556 c Distance restraints
559 C Copy the coordinates from reference coordinates (?)
563 c write (iout,*) "c(",j,i,") =",c(j,i)
567 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
569 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
570 open (ientin,file=tpl_k_rescore,status='old')
571 if (nnt.gt.1) rescore(k,1)=0.0d0
572 do irec=nnt,nct ! loop for reading res sim
574 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
575 & rescore3_tmp,idomain_tmp
577 idomain(k,i_tmp)=idomain_tmp
578 rescore(k,i_tmp)=rescore_tmp
579 rescore2(k,i_tmp)=rescore2_tmp
580 rescore3(k,i_tmp)=rescore3_tmp
581 write(iout,'(a7,i5,3f10.5,i5)') "rescore",
582 & i_tmp,rescore2_tmp,rescore_tmp,
583 & rescore3_tmp,idomain_tmp
586 read (ientin,*,end=1401) rescore_tmp
588 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
589 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
590 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
595 if (waga_dist.ne.0.0d0) then
603 distal=dsqrt(x12*x12+y12*y12+z12*z12)
604 c write (iout,*) k,i,j,distal,dist2_cut
606 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
607 & .and. distal.le.dist2_cut ) then
613 c write (iout,*) "k",k
614 c write (iout,*) "i",i," j",j," constr_homology",
622 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
624 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
625 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
626 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
628 if (odl(k,ii).le.dist_cut) then
629 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
632 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
633 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
635 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
636 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
640 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
650 c Theta, dihedral and SC retraints
652 if (waga_angle.gt.0.0d0) then
653 c open (ientin,file=tpl_k_sigma_dih,status='old')
654 c do irec=1,maxres-3 ! loop for reading sigma_dih
655 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
656 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
657 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
658 c & sigma_dih(k,i+nnt-1)
663 if (idomain(k,i).eq.0) then
667 dih(k,i)=phiref(i) ! right?
668 c read (ientin,*) sigma_dih(k,i) ! original variant
669 c write (iout,*) "dih(",k,i,") =",dih(k,i)
670 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
671 c & "rescore(",k,i-1,") =",rescore(k,i-1),
672 c & "rescore(",k,i-2,") =",rescore(k,i-2),
673 c & "rescore(",k,i-3,") =",rescore(k,i-3)
675 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
676 & rescore(k,i-2)+rescore(k,i-3))/4.0
677 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
678 c write (iout,*) "Raw sigmas for dihedral angle restraints"
679 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
680 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
681 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
682 c Instead of res sim other local measure of b/b str reliability possible
683 if (sigma_dih(k,i).ne.0)
684 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
685 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
690 if (waga_theta.gt.0.0d0) then
691 c open (ientin,file=tpl_k_sigma_theta,status='old')
692 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
693 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
694 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
695 c & sigma_theta(k,i+nnt-1)
700 do i = nnt+2,nct ! right? without parallel.
701 c do i = i=1,nres ! alternative for bounds acc to readpdb?
702 c do i=ithet_start,ithet_end ! with FG parallel.
703 if (idomain(k,i).eq.0) then
707 thetatpl(k,i)=thetaref(i)
708 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
709 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
710 c & "rescore(",k,i-1,") =",rescore(k,i-1),
711 c & "rescore(",k,i-2,") =",rescore(k,i-2)
712 c read (ientin,*) sigma_theta(k,i) ! 1st variant
713 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
714 & rescore(k,i-2))/3.0
715 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
716 if (sigma_theta(k,i).ne.0)
717 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
719 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
720 c rescore(k,i-2) ! right expression ?
721 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
725 if (waga_d.gt.0.0d0) then
726 c open (ientin,file=tpl_k_sigma_d,status='old')
727 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
728 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
729 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
730 c & sigma_d(k,i+nnt-1)
734 do i = nnt,nct ! right? without parallel.
735 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
736 c do i=loc_start,loc_end ! with FG parallel.
737 if (itype(i).eq.10) cycle
738 if (idomain(k,i).eq.0 ) then
745 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
746 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
747 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
748 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
749 sigma_d(k,i)=rescore3(k,i) ! right expression ?
750 if (sigma_d(k,i).ne.0)
751 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
753 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
754 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
755 c read (ientin,*) sigma_d(k,i) ! 1st variant
760 c remove distance restraints not used in any model from the list
761 c shift data in all arrays
763 if (waga_dist.ne.0.0d0) then
769 if (ii_in_use(ii).eq.0.and.liiflag) then
773 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
774 & .not.liiflag.and.ii.eq.lim_odl) then
775 if (ii.eq.lim_odl) then
781 do ki=iistart,lim_odl-iishift
782 ires_homo(ki)=ires_homo(ki+iishift)
783 jres_homo(ki)=jres_homo(ki+iishift)
784 ii_in_use(ki)=ii_in_use(ki+iishift)
785 do k=1,constr_homology
786 odl(k,ki)=odl(k,ki+iishift)
787 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
788 l_homo(k,ki)=l_homo(k,ki+iishift)
792 lim_odl=lim_odl-iishift
798 endif ! .not. klapaucjusz
800 if (constr_homology.gt.0) call homology_partition
801 if (constr_homology.gt.0) call init_int_table
802 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
803 cd & "lim_xx=",lim_xx
804 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
805 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
809 if (.not.lprn) return
810 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
811 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
812 write (iout,*) "Distance restraints from templates"
814 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
815 & ii,ires_homo(ii),jres_homo(ii),
816 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
817 & ki=1,constr_homology)
819 write (iout,*) "Dihedral angle restraints from templates"
821 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
822 & (rad2deg*dih(ki,i),
823 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
825 write (iout,*) "Virtual-bond angle restraints from templates"
827 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
828 & (rad2deg*thetatpl(ki,i),
829 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
831 write (iout,*) "SC restraints from templates"
833 write(iout,'(i5,100(4f8.2,4x))') i,
834 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
835 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
838 c -----------------------------------------------------------------
841 c----------------------------------------------------------------------
842 subroutine read_klapaucjusz
845 include 'DIMENSIONS.ZSCOPT'
846 include 'DIMENSIONS.FREE'
850 include 'COMMON.SETUP'
851 include 'COMMON.CONTROL'
852 include 'COMMON.CHAIN'
853 include 'COMMON.IOUNITS'
855 include 'COMMON.INTERACT'
856 include 'COMMON.NAMES'
857 include 'COMMON.HOMRESTR'
858 character*256 fragfile
859 integer ninclust(maxclust),inclust(max_template,maxclust),
860 & nresclust(maxclust),iresclust(maxres,maxclust)
863 character*24 model_ki_dist, model_ki_angle
864 character*500 controlcard
865 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
866 integer idomain(max_template,maxres)
867 logical lprn /.true./
870 logical unres_pdb,liiflag
873 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
874 double precision, dimension (max_template,maxres) :: rescore
875 double precision, dimension (max_template,maxres) :: rescore2
876 character*24 tpl_k_rescore
883 double precision chomo(3,maxres2+2,max_template)
884 call getenv("FRAGFILE",fragfile)
885 open(ientin,file=fragfile,status="old",err=10)
886 read(ientin,*) constr_homology,nclust
892 do k=1,constr_homology
893 read(ientin,'(a)') pdbfile
894 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
895 & pdbfile(:ilen(pdbfile))
896 open(ipdbin,file=pdbfile,status='old',err=33)
898 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
899 & pdbfile(:ilen(pdbfile))
903 call readpdb_template(k)
916 read(ientin,*) ninclust(i),nresclust(i)
917 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
918 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
924 do ll = 1,ninclust(l)
932 idomain(k,iresclust(i,l)+1) = 1
934 idomain(k,iresclust(i,l)) = 1
938 c Distance restraints
941 C Copy the coordinates from reference coordinates (?)
945 c write (iout,*) "c(",j,i,") =",c(j,i)
948 call int_from_cart1(.false.)
949 call int_from_cart(.true.,.false.)
950 call sc_loc_geom(.false.)
955 if (waga_dist.ne.0.0d0) then
963 distal=dsqrt(x12*x12+y12*y12+z12*z12)
964 c write (iout,*) k,i,j,distal,dist2_cut
966 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
967 & .and. distal.le.dist2_cut ) then
973 c write (iout,*) "k",k
974 c write (iout,*) "i",i," j",j," constr_homology",
982 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
984 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
985 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
986 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
988 if (odl(k,ii).le.dist_cut) then
989 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
992 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
993 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
995 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
996 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1000 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1003 c l_homo(k,ii)=.false.
1010 c Theta, dihedral and SC retraints
1012 if (waga_angle.gt.0.0d0) then
1014 if (idomain(k,i).eq.0) then
1015 c sigma_dih(k,i)=0.0
1019 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
1020 & rescore(k,i-2)+rescore(k,i-3))/4.0
1021 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
1022 c & " sigma_dihed",sigma_dih(k,i)
1023 if (sigma_dih(k,i).ne.0)
1024 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1029 if (waga_theta.gt.0.0d0) then
1031 if (idomain(k,i).eq.0) then
1032 c sigma_theta(k,i)=0.0
1035 thetatpl(k,i)=thetaref(i)
1036 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
1037 & rescore(k,i-2))/3.0
1038 if (sigma_theta(k,i).ne.0)
1039 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1043 if (waga_d.gt.0.0d0) then
1045 if (itype(i).eq.10) cycle
1046 if (idomain(k,i).eq.0 ) then
1053 sigma_d(k,i)=rescore(k,i)
1054 if (sigma_d(k,i).ne.0)
1055 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1056 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
1062 c remove distance restraints not used in any model from the list
1063 c shift data in all arrays
1065 if (waga_dist.ne.0.0d0) then
1071 if (ii_in_use(ii).eq.0.and.liiflag) then
1075 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
1076 & .not.liiflag.and.ii.eq.lim_odl) then
1077 if (ii.eq.lim_odl) then
1078 iishift=ii-iistart+1
1083 do ki=iistart,lim_odl-iishift
1084 ires_homo(ki)=ires_homo(ki+iishift)
1085 jres_homo(ki)=jres_homo(ki+iishift)
1086 ii_in_use(ki)=ii_in_use(ki+iishift)
1087 do k=1,constr_homology
1088 odl(k,ki)=odl(k,ki+iishift)
1089 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
1090 l_homo(k,ki)=l_homo(k,ki+iishift)
1094 lim_odl=lim_odl-iishift
1101 10 stop "Error infragment file"
1103 c----------------------------------------------------------------------