subroutine read_control C C Read molecular data C implicit none include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.HEADER' include 'COMMON.FFIELD' include 'COMMON.INTERACT' include "COMMON.SPLITELE" include 'COMMON.SHIELD' include 'COMMON.FREE' include 'COMMON.LANGEVIN' character*320 controlcard,ucase #ifdef MPL include 'COMMON.INFO' #endif integer i,i1,i2,it1,it2 double precision pi read (INP,'(a80)') titel call card_concat(controlcard) call readi(controlcard,'NRES',nres,0) call readi(controlcard,'RESCALE',rescale_mode,2) call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) write (iout,*) "DISTCHAINMAX",distchainmax C Reading the dimensions of box in x,y,z coordinates call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) call reada(controlcard,"LIPTHICK",lipthick,0.0d0) call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) if (lipthick.gt.0.0d0) then bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick C endif if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick if ((lipbufthick*2.0d0).gt.lipthick) &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" endif write(iout,*) "bordliptop=",bordliptop write(iout,*) "bordlipbot=",bordlipbot write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot C Shielding mode call readi(controlcard,'SHIELD',shield_mode,0) write (iout,*) "SHIELD MODE",shield_mode if (shield_mode.gt.0) then pi=3.141592d0 C VSolvSphere the volume of solving sphere C print *,pi,"pi" C rpp(1,1) is the energy r0 for peptide group contact and will be used for it C there will be no distinction between proline peptide group and normal peptide C group in case of shielding parameters VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 write (iout,*) VSolvSphere,VSolvSphere_div C long axis of side chain do i=1,ntyp long_r_sidechain(i)=vbldsc0(1,i) short_r_sidechain(i)=sigma0(i) enddo buff_shield=1.0d0 endif call readi(controlcard,'PDBOUT',outpdb,0) call readi(controlcard,'MOL2OUT',outmol2,0) refstr=(index(controlcard,'REFSTR').gt.0) write (iout,*) "REFSTR",refstr pdbref=(index(controlcard,'PDBREF').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 call readi(controlcard,'CONSTR_DIST',constr_dist,0) write (iout,*) "with_dihed_constr ",with_dihed_constr, & " CONSTR_DIST",constr_dist with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 write (iout,*) "with_theta_constr ",with_theta_constr call flush(iout) min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) call readi(controlcard,'NCUT',ncut,0) call readi(controlcard,'NCLUST',nclust,5) call readi(controlcard,'SYM',symetr,1) write (iout,*) 'sym', symetr call readi(controlcard,'NSTART',nstart,0) call readi(controlcard,'NEND',nend,0) call reada(controlcard,'ECUT',ecut,10.0d0) call reada(controlcard,'PROB',prob_limit,0.99d0) write (iout,*) "Probability limit",prob_limit lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) if (ncut.gt.0) & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) call readi(controlcard,'IOPT',iopt,2) lside = index(controlcard,"SIDE").gt.0 efree = index(controlcard,"EFREE").gt.0 call readi(controlcard,'NTEMP',nT,1) write (iout,*) "nT",nT call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0) write (iout,*) "nT",nT write (iout,*) 'beta_h',(beta_h(i),i=1,nT) do i=1,nT beta_h(i)=1.0d0/(1.987D-3*beta_h(i)) enddo write (iout,*) 'beta_h',(beta_h(i),i=1,nT) lprint_cart=index(controlcard,"PRINT_CART") .gt.0 lprint_int=index(controlcard,"PRINT_INT") .gt.0 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) write (iout,*) "with_homology_constr ",with_dihed_constr, & " CONSTR_HOMOLOGY",constr_homology print_homology_restraints= & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0 print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0 print_homology_models= & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0 call readi(controlcard,'NSAXS',nsaxs,0) call readi(controlcard,'SAXS_MODE',saxs_mode,0) call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0) call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0) write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE", & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff if (min_var) iopt=1 return end c-------------------------------------------------------------------------- subroutine molread C C Read molecular data. C implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.CONTACTS' include 'COMMON.TIME1' include 'COMMON.TORCNSTR' include 'COMMON.SHIELD' #ifdef MPL include 'COMMON.INFO' #endif character*4 sequence(maxres) character*800 weightcard integer rescode double precision x(maxvar) integer itype_pdb(maxres) logical seq_comp integer i,j,kkk,i1,i2,it1,it2 C C Body C C Read weights of the subsequent energy terms. call card_concat(weightcard) write(iout,*) weightcard call reada(weightcard,'WSC',wsc,1.0d0) call reada(weightcard,'WLONG',wsc,wsc) call reada(weightcard,'WSCP',wscp,1.0d0) call reada(weightcard,'WELEC',welec,1.0D0) call reada(weightcard,'WVDWPP',wvdwpp,welec) call reada(weightcard,'WEL_LOC',wel_loc,1.0D0) call reada(weightcard,'WCORR4',wcorr4,0.0D0) call reada(weightcard,'WCORR5',wcorr5,0.0D0) call reada(weightcard,'WCORR6',wcorr6,0.0D0) call reada(weightcard,'WTURN3',wturn3,1.0D0) call reada(weightcard,'WTURN4',wturn4,1.0D0) call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'WSAXS',wsaxs,0.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) if (index(weightcard,'SOFT').gt.0) ipot=6 call reada(weightcard,"D0CM",d0cm,3.78d0) call reada(weightcard,"AKCM",akcm,15.1d0) call reada(weightcard,"AKTH",akth,11.0d0) call reada(weightcard,"AKCT",akct,12.0d0) call reada(weightcard,"V1SS",v1ss,-1.08d0) call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) call reada(weightcard,'WSHIELD',wshield,1.0d0) write(iout,*) 'WSHIELD',wshield call reada(weightcard,'WLT',wliptran,0.0D0) call reada(weightcard,"ATRISS",atriss,0.301D0) call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) call reada(weightcard,"DTRISS",dtriss,1.001D0) write (iout,*) "ATRISS=", atriss write (iout,*) "BTRISS=", btriss write (iout,*) "CTRISS=", ctriss write (iout,*) "DTRISS=", dtriss dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. enddo do i=1,maxres-1 do j=i+1,maxres dyn_ssbond_ij(i,j)=1.0d300 enddo enddo call reada(weightcard,"HT",Ht,0.0D0) if (dyn_ss) then ss_depth=ebr/wsc-0.25*eps(1,1) Ht=Ht/wsc-0.25*eps(1,1) akcm=akcm*wstrain/wsc akth=akth*wstrain/wsc akct=akct*wstrain/wsc v1ss=v1ss*wstrain/wsc v2ss=v2ss*wstrain/wsc v3ss=v3ss*wstrain/wsc else ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain endif write (iout,'(/a)') "Disulfide bridge parameters:" write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, & ' v3ss:',v3ss C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc weights(2)=wscp weights(3)=welec weights(4)=wcorr weights(5)=wcorr5 weights(6)=wcorr6 weights(7)=wel_loc weights(8)=wturn3 weights(9)=wturn4 weights(10)=wturn6 weights(11)=wang weights(12)=wscloc weights(13)=wtor weights(14)=wtor_d weights(15)=wstrain weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 weights(19)=wsccor write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3, & wturn4,wturn6,wsccor 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ & 'WELEC= ',f10.6,' (p-p electr)'/ & 'WVDWPP= ',f10.6,' (p-p VDW)'/ & 'WBOND= ',f10.6,' (stretching)'/ & 'WANG= ',f10.6,' (bending)'/ & 'WSCLOC= ',f10.6,' (SC local)'/ & 'WTOR= ',f10.6,' (torsional)'/ & 'WTORD= ',f10.6,' (double torsional)'/ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/ & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)'/ & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)') if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' write (iout,'(2(a,f5.3/))') & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr, & 'Range of quenching the correlation terms:',2*delt_corr else if (wcorr.gt.0.0d0) then write (iout,'(/2a/)') 'Hydrogen-bonding correlation ', & 'between contact pairs of peptide groups' endif write (iout,'(a,f8.3)') & 'Scaling factor of 1,4 SC-p interactions:',scal14 write (iout,'(a,f8.3)') & 'General scaling factor of SC-p interactions:',scalscp r0_corr=cutoff_corr-delt_corr do i=1,20 aad(i,1)=scalscp*aad(i,1) aad(i,2)=scalscp*aad(i,2) bad(i,1)=scalscp*bad(i,1) bad(i,2)=scalscp*bad(i,2) enddo call flush(iout) print *,'indpdb=',indpdb,' pdbref=',pdbref C Read sequence if not taken from the pdb file. if (iscode.gt.0) then read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) else read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) endif C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo print *,nres print '(20i4)',(itype(i),i=1,nres) do i=1,nres #ifdef PROCOR if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR else if (iabs(itype(i+1)).ne.20) then #else else if (iabs(itype(i)).ne.20) then #endif itel(i)=1 else itel(i)=2 endif enddo write (iout,*) "ITEL" do i=1,nres-1 write (iout,*) i,itype(i),itel(i) enddo print *,'Call Read_Bridge.' call read_bridge C this fragment reads diheadral constrains if (with_dihed_constr) then read (inp,*) ndih_constr if (ndih_constr.gt.0) then C read (inp,*) ftors C write (iout,*) 'FTORS',ftors C ftors is the force constant for torsional quartic constrains read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), & i=1,ndih_constr) write (iout,*) & 'There are',ndih_constr,' constraints on phi angles.' do i=1,ndih_constr write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), & ftors(i) enddo do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo endif ! endif ndif_constr.gt.0 endif ! with_dihed_constr if (with_theta_constr) then C with_theta_constr is keyword allowing for occurance of theta constrains read (inp,*) ntheta_constr C ntheta_constr is the number of theta constrains if (ntheta_constr.gt.0) then C read (inp,*) ftors read (inp,*) (itheta_constr(i),theta_constr0(i), & theta_drange(i),for_thet_constr(i), & i=1,ntheta_constr) C the above code reads from 1 to ntheta_constr C itheta_constr(i) residue i for which is theta_constr C theta_constr0 the global minimum value C theta_drange is range for which there is no energy penalty C for_thet_constr is the force constant for quartic energy penalty C E=k*x**4 C if(me.eq.king.or..not.out1file)then write (iout,*) & 'There are',ntheta_constr,' constraints on phi angles.' do i=1,ntheta_constr write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), & theta_drange(i), & for_thet_constr(i) enddo C endif do i=1,ntheta_constr theta_constr0(i)=deg2rad*theta_constr0(i) theta_drange(i)=deg2rad*theta_drange(i) enddo C if(me.eq.king.or..not.out1file) C & write (iout,*) 'FTORS',ftors C do i=1,ntheta_constr C ii = itheta_constr(i) C thetabound(1,ii) = phi0(i)-drange(i) C thetabound(2,ii) = phi0(i)+drange(i) C enddo endif ! ntheta_constr.gt.0 endif! with_theta_constr nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.ntyp1) nnt=2 if (itype(nres).eq.ntyp1) nct=nct-1 if (nstart.lt.nnt) nstart=nnt if (nend.gt.nct .or. nend.eq.0) nend=nct write (iout,*) "nstart",nstart," nend",nend write (iout,*) "calling read_saxs_consrtr",nsaxs if (nsaxs.gt.0) call read_saxs_constr nres0=nres if (constr_homology.gt.0) then call read_constr_homology(print_homology_restraints) endif c if (pdbref) then c read(inp,'(a)') pdbfile c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile c open(ipdbin,file=pdbfile,status='old',err=33) c goto 34 c 33 write (iout,'(a)') 'Error opening PDB file.' c stop c 34 continue c print *,'Begin reading pdb data' c call readpdb c print *,'Finished reading pdb data' c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup c do i=1,nres c itype_pdb(i)=itype(i) c enddo c close (ipdbin) c write (iout,'(a,i3)') 'nsup=',nsup c nstart_seq=nnt c if (nsup.le.(nct-nnt+1)) then c do i=0,nct-nnt+1-nsup c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then c nstart_seq=nnt+i c goto 111 c endif c enddo c write (iout,'(a)') c & 'Error - sequences to be superposed do not match.' c stop c else c do i=0,nsup-(nct-nnt+1) c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) c & then c nstart_sup=nstart_sup+i c nsup=nct-nnt+1 c goto 111 c endif c enddo c write (iout,'(a)') c & 'Error - sequences to be superposed do not match.' c endif c 111 continue c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup, c & ' nstart_seq=',nstart_seq c endif call init_int_table call setup_var write (iout,*) "molread: REFSTR",refstr if (refstr) then if (.not.pdbref) then call read_angles(inp,*38) goto 39 38 write (iout,'(a)') 'Error reading reference structure.' #ifdef MPL call mp_stopall(Error_Msg) #else stop 'Error reading reference structure' #endif 39 call chainbuild nstart_sup=nnt nstart_seq=nnt nsup=nct-nnt+1 kkk=1 do i=1,2*nres do j=1,3 cref(j,i,kkk)=c(j,i) enddo enddo endif call contact(.true.,ncont_ref,icont_ref) endif if (ns.gt.0) then C write (iout,'(/a,i3,a)') C & 'The chain contains',ns,' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) if (dyn_ss) then write(iout,*)"Running with dynamic disulfide-bond formation" else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres i2=jhpb(i)-nres it1=itype(i1) it2=itype(i2) write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i), & ebr,forcon(i) enddo write (iout,'(a)') endif endif if (ns.gt.0.and.dyn_ss) then do i=nss+1,nhpb ihpb(i-nss)=ihpb(i) jhpb(i-nss)=jhpb(i) forcon(i-nss)=forcon(i) dhpb(i-nss)=dhpb(i) enddo nhpb=nhpb-nss nss=0 call hpb_partition do i=1,ns dyn_ss_mask(iss(i))=.true. enddo endif c Read distance restraints if (constr_dist.gt.0) then call read_dist_constr call hpb_partition endif return end c----------------------------------------------------------------------------- logical function seq_comp(itypea,itypeb,length) implicit none integer length,itypea(length),itypeb(length) integer i do i=1,length if (itypea(i).ne.itypeb(i)) then seq_comp=.false. return endif enddo seq_comp=.true. return end c----------------------------------------------------------------------------- subroutine read_bridge C Read information about disulfide bridges. implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' include 'COMMON.TIME1' #ifdef MPL include 'COMMON.INFO' #endif integer i,j C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) print *,'ns=',ns C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then write (iout,'(2a,i3,a)') & 'Do you REALLY think that the residue ', & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ', & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' #ifdef MPL call mp_stopall(error_msg) #else stop #endif endif enddo C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of C bridging residues. do i=1,nss do j=1,i-1 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then write (iout,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' write (*,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' #ifdef MPL call mp_stopall(error_msg) #else stop #endif endif enddo do j=1,ns if (ihpb(i).eq.iss(j)) goto 10 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 10 continue do j=1,ns if (jhpb(i).eq.iss(j)) goto 20 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 20 continue C dhpb(i)=dbr C forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres enddo endif endif return end c---------------------------------------------------------------------------- subroutine read_angles(kanal,*) implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' integer i,kanal read (kanal,*,err=10,end=10) (theta(i),i=3,nres) read (kanal,*,err=10,end=10) (phi(i),i=4,nres) read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1) read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1) do i=1,nres theta(i)=deg2rad*theta(i) phi(i)=deg2rad*phi(i) alph(i)=deg2rad*alph(i) omeg(i)=deg2rad*omeg(i) enddo return 10 return1 end c---------------------------------------------------------------------------- subroutine reada(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch double precision wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end c---------------------------------------------------------------------------- subroutine multreada(rekord,lancuch,tablica,dim,default) implicit none integer dim,i double precision tablica(dim),default character*(*) rekord,lancuch integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,lancuch) if (iread.eq.0) return iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end c---------------------------------------------------------------------------- subroutine readi(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch integer wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end C---------------------------------------------------------------------- subroutine multreadi(rekord,lancuch,tablica,dim,default) implicit none integer dim,i integer tablica(dim),default character*(*) rekord,lancuch character*80 aux integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) return iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end c---------------------------------------------------------------------------- subroutine card_concat(card) include 'DIMENSIONS' include 'COMMON.IOUNITS' character*(*) card character*80 karta,ucase external ilen read (inp,'(a)') karta karta=ucase(karta) card=' ' do while (karta(80:80).eq.'&') card=card(:ilen(card)+1)//karta(:79) read (inp,'(a)') karta karta=ucase(karta) enddo card=card(:ilen(card)+1)//karta return end c---------------------------------------------------------------------------- subroutine openunits implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" character*3 liczba include "COMMON.MPI" #endif include 'COMMON.IOUNITS' include 'COMMON.CONTROL' integer lenpre,lenpot,ilen external ilen character*16 cformat,cprint character*16 ucase integer lenint,lenout call getenv('INPUT',prefix) call getenv('OUTPUT',prefout) call getenv('INTIN',prefintin) call getenv('COORD',cformat) call getenv('PRINTCOOR',cprint) call getenv('SCRATCHDIR',scratchdir) from_bx=.true. from_cx=.false. if (index(ucase(cformat),'CX').gt.0) then from_cx=.true. from_bx=.false. endif from_cart=.true. lenpre=ilen(prefix) lenout=ilen(prefout) lenint=ilen(prefintin) C Get the names and open the input files open (inp,file=prefix(:ilen(prefix))//'.inp',status='old') #ifdef MPI write (liczba,'(bz,i3.3)') me outname=prefout(:lenout)//'_clust.out_'//liczba #else outname=prefout(:lenout)//'_clust.out' #endif if (from_bx) then intinname=prefintin(:lenint)//'.bx' else if (from_cx) then intinname=prefintin(:lenint)//'.cx' else intinname=prefintin(:lenint)//'.int' endif rmsname=prefintin(:lenint)//'.rms' open (jplot,file=prefout(:ilen(prefout))//'.tex', & status='unknown') open (jrms,file=rmsname,status='unknown') open(iout,file=outname,status='unknown') C Get parameter filenames and open the parameter files. call getenv('BONDPAR',bondname) open (ibond,file=bondname,status='old') call getenv('THETPAR',thetname) open (ithep,file=thetname,status='old') call getenv('ROTPAR',rotname) open (irotam,file=rotname,status='old') call getenv('TORPAR',torname) open (itorp,file=torname,status='old') call getenv('TORDPAR',tordname) open (itordp,file=tordname,status='old') call getenv('FOURIER',fouriername) open (ifourier,file=fouriername,status='old') call getenv('ELEPAR',elename) open (ielep,file=elename,status='old') call getenv('SIDEPAR',sidename) open (isidep,file=sidename,status='old') call getenv('SIDEP',sidepname) open (isidep1,file=sidepname,status="old") call getenv('SCCORPAR',sccorname) open (isccor,file=sccorname,status="old") call getenv('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old') #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file C Use -DOLDSCP to use hard-coded constants instead. C call getenv('SCPPAR',scpname) open (iscpp,file=scpname,status='old') #endif return end subroutine read_dist_constr implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.SBRIDGE' integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard logical lprn /.true./ logical normalize,next integer restr_type double precision xlink(4,0:4) / c a b c sigma & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) next=.true. DO WHILE (next) call card_concat(controlcard) next = index(controlcard,"NEXT").gt.0 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist) write (iout,*) "restr_type",restr_type call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) call readi(controlcard,"NDIST",ndist_,0) call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) if (restr_type.eq.10) & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0) call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0) call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) normalize = index(controlcard,"NORMALIZE").gt.0 write (iout,*) "WBOLTZD",wboltzd c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ c write (iout,*) "IFRAG" c do i=1,nfrag_ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) c enddo c write (iout,*) "IPAIR" c do i=1,npair_ c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) c enddo if (nfrag_.gt.0) then nres0=nres read(inp,'(a)') pdbfile write (iout,*) & "Distance restraints will be constructed from structure ",pdbfile open(ipdbin,file=pdbfile,status='old',err=11) call readpdb(.true.) nres=nres0 close(ipdbin) endif do i=1,nfrag_ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup if (ifrag_(2,i).gt.nstart_sup+nsup-1) & ifrag_(2,i)=nstart_sup+nsup-1 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) c call flush(iout) if (wfrag_(i).eq.0.0d0) cycle do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) c write (iout,*) "j",j," k",k ddjk=dist(j,k) if (restr_type.eq.1) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) endif else if (restr_type.eq.3) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) endif write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) enddo enddo enddo do i=1,npair_ if (wpair_(i).eq.0.0d0) cycle ii = ipair_(1,i) jj = ipair_(2,i) if (ii.gt.jj) then itemp=ii ii=jj jj=itemp endif do j=ifrag_(1,ii),ifrag_(2,ii) do k=ifrag_(1,jj),ifrag_(2,jj) if (restr_type.eq.1) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) endif else if (restr_type.eq.3) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) endif write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) enddo enddo enddo c print *,ndist_ write (iout,*) "Distance restraints as read from input" do i=1,ndist_ if (restr_type.eq.11) then read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle nhpb=nhpb+1 irestr_type(nhpb)=11 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) if (ibecarb(nhpb).gt.0) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif else if (constr_dist.eq.10) then c Cross-lonk Markov-like potential call card_concat(controlcard) call readi(controlcard,"ILINK",ihpb(nhpb+1),0) call readi(controlcard,"JLINK",jhpb(nhpb+1),0) ibecarb(nhpb+1)=0 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle if (index(controlcard,"ZL").gt.0) then link_type=1 else if (index(controlcard,"ADH").gt.0) then link_type=2 else if (index(controlcard,"PDH").gt.0) then link_type=3 else if (index(controlcard,"DSS").gt.0) then link_type=4 else link_type=0 endif call reada(controlcard,"AXLINK",dhpb(nhpb+1), & xlink(1,link_type)) call reada(controlcard,"BXLINK",dhpb1(nhpb+1), & xlink(2,link_type)) call reada(controlcard,"CXLINK",fordepth(nhpb+1), & xlink(3,link_type)) call reada(controlcard,"SIGMA",forcon(nhpb+1), & xlink(4,link_type)) call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0) c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1), c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) if (forcon(nhpb+1).le.0.0d0 .or. & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle nhpb=nhpb+1 irestr_type(nhpb)=10 if (ibecarb(nhpb).gt.0) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) else C print *,"in else" read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 if (dhpb1(nhpb).eq.0.0d0) then irestr_type(nhpb)=1 else irestr_type(nhpb)=2 endif if (ibecarb(nhpb).gt.0) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif if (dhpb(nhpb).eq.0.0d0) & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) endif write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) endif C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) C if (forcon(nhpb+1).gt.0.0d0) then C nhpb=nhpb+1 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) enddo ENDDO ! next fordepthmax=0.0d0 if (normalize) then do i=nss+1,nhpb if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax) & fordepthmax=fordepth(i) enddo do i=nss+1,nhpb if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax enddo endif if (nhpb.gt.nss) then write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)') & "The following",nhpb-nss, & " distance restraints have been imposed:", & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V", & " score"," type" do i=nss+1,nhpb write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i), & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i), & irestr_type(i) enddo endif call hpb_partition call flush(iout) return 11 write (iout,*)"read_dist_restr: error reading reference structure" stop end c------------------------------------------------------------------------------- subroutine read_saxs_constr implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.SBRIDGE' double precision cm(3) c read(inp,*) nsaxs write (iout,*) "Calling read_saxs nsaxs",nsaxs call flush(iout) if (saxs_mode.eq.0) then c SAXS distance distribution do i=1,nsaxs read(inp,*) distsaxs(i),Psaxs(i) enddo Cnorm = 0.0d0 do i=1,nsaxs Cnorm = Cnorm + Psaxs(i) enddo write (iout,*) "Cnorm",Cnorm do i=1,nsaxs Psaxs(i)=Psaxs(i)/Cnorm enddo write (iout,*) "Normalized distance distribution from SAXS" do i=1,nsaxs write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i) enddo Wsaxs0=0.0d0 do i=1,nsaxs Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i)) enddo write (iout,*) "Wsaxs0",Wsaxs0 else c SAXS "spheres". do i=1,nsaxs read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3) enddo do j=1,3 cm(j)=0.0d0 enddo do i=1,nsaxs do j=1,3 cm(j)=cm(j)+Csaxs(j,i) enddo enddo do j=1,3 cm(j)=cm(j)/nsaxs enddo do i=1,nsaxs do j=1,3 Csaxs(j,i)=Csaxs(j,i)-cm(j) enddo enddo write (iout,*) "SAXS sphere coordinates" do i=1,nsaxs write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3) enddo endif return end c====------------------------------------------------------------------- subroutine read_constr_homology(lprn) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.HOMRESTR' c c For new homol impl c include 'COMMON.VAR' c include 'include_unres/COMMON.VAR' c c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d, c & dist_cut c common /przechowalnia/ odl_temp(maxres,maxres,max_template), c & sigma_odl_temp(maxres,maxres,max_template) character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp integer idomain(max_template,maxres) integer ilen external ilen logical lprn logical unres_pdb,liiflag c c FP - Nov. 2014 Temporary specifications for new vars c double precision rescore_tmp,x12,y12,z12,rescore2_tmp, & rescore3_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 double precision, dimension (max_template,maxres) :: rescore3 character*24 tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints c not using pre-computed ones stored in files model_ki_{dist,angle} c FP (Nov., 2014) c ----------------------------------------------------------------- c c c Alternative: reading from input #ifdef DEBUG write (iout,*) "BEGIN READ HOMOLOGY INFO" #ifdef AIX call flush_(iout) #else call flush(iout) #endif #endif call card_concat(controlcard) call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0) call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0) call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) dist1cut=(index(controlcard,'DIST1CUT').gt.0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) if (homol_nset.gt.1)then call readi(controlcard,"ISET",iset,1) call card_concat(controlcard) read(controlcard,*) (waga_homology(i),i=1,homol_nset) else iset=1 waga_homology(1)=1.0 endif c #ifdef DEBUG write(iout,*) "read_constr_homology iset",iset write(iout,*) "waga_homology(",iset,")",waga_homology(iset) #ifdef AIX call flush_(iout) #else call flush(iout) #endif #endif cd write (iout,*) "nnt",nnt," nct",nct cd call flush(iout) lim_odl=0 lim_dih=0 c c New c c c Reading HM global scores (prob not required) c do i = nnt,nct do k=1,constr_homology idomain(k,i)=0 enddo enddo c open (4,file="HMscore") c do k=1,constr_homology c read (4,*,end=521) hmscore_tmp c hmscore(k)=hmscore_tmp ! Another transformation can be used c write(*,*) "Model", k, ":", hmscore(k) c enddo c521 continue ii=0 do i = nnt,nct-2 do j=i+2,nct ii=ii+1 ii_in_use(ii)=0 enddo enddo c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d if (read_homol_frag) then call read_klapaucjusz else write (iout,*) "CONSTR_HOMOLOGY",constr_homology do k=1,constr_homology read(inp,'(a)') pdbfile c write (iout,*) "k ",k," pdbfile ",pdbfile c Next stament causes error upon compilation (?) c if(me.eq.king.or. .not. out1file) c write (iout,'(2a)') 'PDB data will be read from file ', c & pdbfile(:ilen(pdbfile)) write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file', & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file', & pdbfile(:ilen(pdbfile)) stop 34 continue c print *,'Begin reading pdb data' c c Files containing res sim or local scores (former containing sigmas) c write(kic2,'(bz,i2.2)') k tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. call readpdb_template(k) do i=1,2*nres do j=1,3 crefjlee(j,i)=c(j,i) enddo enddo #ifdef DEBUG do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo write (iout,*) "READ HOMOLOGY INFO" write (iout,*) "read_constr_homology x: after reading pdb file" write (iout,*) "waga_homology(",iset,")",waga_homology(iset) write (iout,*) "waga_dist",waga_dist write (iout,*) "waga_angle",waga_angle write (iout,*) "waga_theta",waga_theta write (iout,*) "waga_d",waga_d write (iout,*) "dist_cut",dist_cut #endif #ifdef AIX call flush_(iout) #else call flush(iout) #endif c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) do i=1,2*nres do j=1,3 c c(j,i)=cref(j,i) c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo c c From read_dist_constr (commented out 25/11/2014 <-> res sim) c c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore open (ientin,file=tpl_k_rescore,status='old') if (nnt.gt.1) rescore(k,1)=0.0d0 do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, & rescore3_tmp,idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp rescore3(k,i_tmp)=rescore3_tmp write(iout,'(a7,i5,3f10.5,i5)') "rescore", & i_tmp,rescore2_tmp,rescore_tmp, & rescore3_tmp,idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) endif enddo 1401 continue close (ientin) if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 do j=i+2,nct x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut if (dist1cut .and. k.gt.1) then ii=ii+1 if (l_homo(1,ii)) then ii_in_use(ii)=1 l_homo(k,ii)=.true. ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal sigma_odl(k,ii)=sigma_odl(1,ii) else l_homo(k,ii)=.false. endif else if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", c & constr_homology ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 l_homo(k,ii)=.false. endif endif enddo enddo lim_odl=ii endif c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_dih,status='old') c do irec=1,maxres-3 ! loop for reading sigma_dih c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for? c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right? c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_dih(k,i+nnt-1) c enddo c1402 continue c close (ientin) do i = nnt+3,nct if (idomain(k,i).eq.0) then sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) ! right? c read (ientin,*) sigma_dih(k,i) ! original variant c write (iout,*) "dih(",k,i,") =",dih(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2), c & "rescore(",k,i-3,") =",rescore(k,i-3) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 c write (iout,*) "Raw sigmas for dihedral angle restraints" c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2)*rescore(k,i-3) ! right expression ? c Instead of res sim other local measure of b/b str reliability possible if (sigma_dih(k,i).ne.0) & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) enddo lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_theta,status='old') c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for? c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_theta(k,i+nnt-1) c enddo c1403 continue c close (ientin) do i = nnt+2,nct ! right? without parallel. c do i = i=1,nres ! alternative for bounds acc to readpdb? c do i=ithet_start,ithet_end ! with FG parallel. if (idomain(k,i).eq.0) then sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2) c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 if (sigma_theta(k,i).ne.0) & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for? c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_d(k,i+nnt-1) c enddo c1404 continue do i = nnt,nct ! right? without parallel. c do i=2,nres-1 ! alternative for bounds acc to readpdb? c do i=loc_start,loc_end ! with FG parallel. if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i) sigma_d(k,i)=rescore3(k,i) ! right expression ? if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) c read (ientin,*) sigma_d(k,i) ! 1st variant enddo endif enddo c c remove distance restraints not used in any model from the list c shift data in all arrays c if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 if (ii_in_use(ii).eq.0.and.liiflag) then liiflag=.false. iistart=ii endif if (ii_in_use(ii).ne.0.and..not.liiflag.or. & .not.liiflag.and.ii.eq.lim_odl) then if (ii.eq.lim_odl) then iishift=ii-iistart+1 else iishift=ii-iistart endif liiflag=.true. do ki=iistart,lim_odl-iishift ires_homo(ki)=ires_homo(ki+iishift) jres_homo(ki)=jres_homo(ki+iishift) ii_in_use(ki)=ii_in_use(ki+iishift) do k=1,constr_homology odl(k,ki)=odl(k,ki+iishift) sigma_odl(k,ki)=sigma_odl(k,ki+iishift) l_homo(k,ki)=l_homo(k,ki+iishift) enddo enddo ii=ii-iishift lim_odl=lim_odl-iishift endif enddo enddo endif endif ! .not. klapaucjusz if (constr_homology.gt.0) call homology_partition if (constr_homology.gt.0) call init_int_table cd write (iout,*) "homology_partition: lim_theta= ",lim_theta, cd & "lim_xx=",lim_xx c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end c c Print restraints c if (.not.lprn) return cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" do ii=1,lim_odl write(iout,'(3i5,100(2f8.2,1x,l1,4x))') & ii,ires_homo(ii),jres_homo(ii), & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii), & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" do i=nnt+3,nct write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*dih(ki,i), & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" do i=nnt+2,nct write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*thetatpl(ki,i), & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) enddo endif c ----------------------------------------------------------------- return end c---------------------------------------------------------------------- subroutine read_klapaucjusz include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.SETUP' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.HOMRESTR' character*256 fragfile integer ninclust(maxclust),inclust(max_template,maxclust), & nresclust(maxclust),iresclust(maxres,maxclust) character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp integer idomain(max_template,maxres) logical lprn /.true./ integer ilen external ilen logical unres_pdb,liiflag c c double precision rescore_tmp,x12,y12,z12,rescore2_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 character*24 tpl_k_rescore c c For new homol impl c include 'COMMON.VAR' c double precision chomo(3,maxres2+2,max_template) call getenv("FRAGFILE",fragfile) open(ientin,file=fragfile,status="old",err=10) read(ientin,*) constr_homology,nclust l_homo = .false. sigma_theta=0.0 sigma_d=0.0 sigma_dih=0.0 c Read pdb files do k=1,constr_homology read(ientin,'(a)') pdbfile write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', & pdbfile(:ilen(pdbfile)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file', & pdbfile(:ilen(pdbfile)) stop 34 continue unres_pdb=.false. call readpdb_template(k) do i=1,2*nres do j=1,3 chomo(j,i,k)=c(j,i) enddo enddo do i=1,nres rescore(k,i)=0.2d0 rescore2(k,i)=1.0d0 enddo enddo c Read clusters do i=1,nclust read(ientin,*) ninclust(i),nresclust(i) read(ientin,*) (inclust(k,i),k=1,ninclust(i)) read(ientin,*) (iresclust(k,i),k=1,nresclust(i)) enddo c c Loop over clusters c do l=1,nclust do ll = 1,ninclust(l) k = inclust(ll,l) do i=1,nres idomain(k,i)=0 enddo do i=1,nresclust(l) if (nnt.gt.1) then idomain(k,iresclust(i,l)+1) = 1 else idomain(k,iresclust(i,l)) = 1 endif enddo c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo call int_from_cart1(.false.) call int_from_cart(.true.,.false.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) enddo if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 do j=i+2,nct x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", c & constr_homology ires_homo(ii)=i jres_homo(ii)=j odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 do ik=i,j sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) enddo sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 c l_homo(k,ii)=.false. endif enddo enddo lim_odl=ii endif c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then do i = nnt+3,nct if (idomain(k,i).eq.0) then c sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), c & " sigma_dihed",sigma_dih(k,i) if (sigma_dih(k,i).ne.0) & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) enddo lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then do i = nnt+2,nct if (idomain(k,i).eq.0) then c sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 if (sigma_theta(k,i).ne.0) & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) enddo endif if (waga_d.gt.0.0d0) then do i = nnt,nct if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then c sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) sigma_d(k,i)=rescore(k,i) if (sigma_d(k,i).ne.0) & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 enddo endif enddo ! l enddo ! ll c c remove distance restraints not used in any model from the list c shift data in all arrays c if (waga_dist.ne.0.0d0) then ii=0 liiflag=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 if (ii_in_use(ii).eq.0.and.liiflag) then liiflag=.false. iistart=ii endif if (ii_in_use(ii).ne.0.and..not.liiflag.or. & .not.liiflag.and.ii.eq.lim_odl) then if (ii.eq.lim_odl) then iishift=ii-iistart+1 else iishift=ii-iistart endif liiflag=.true. do ki=iistart,lim_odl-iishift ires_homo(ki)=ires_homo(ki+iishift) jres_homo(ki)=jres_homo(ki+iishift) ii_in_use(ki)=ii_in_use(ki+iishift) do k=1,constr_homology odl(k,ki)=odl(k,ki+iishift) sigma_odl(k,ki)=sigma_odl(k,ki+iishift) l_homo(k,ki)=l_homo(k,ki+iishift) enddo enddo ii=ii-iishift lim_odl=lim_odl-iishift endif enddo enddo endif return 10 stop "Error infragment file" end c----------------------------------------------------------------------