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.FREE' include 'COMMON.INTERACT' include "COMMON.SPLITELE" include 'COMMON.SHIELD' include 'COMMON.SAXS' 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) energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call readi(controlcard,'TORMODE',tor_mode,0) 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,25.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) write (iout,*) "Cutoff on interactions",r_cut write (iout,*) "lambda",rlamb 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 pdbref=(index(controlcard,'PDBREF').gt.0) if (index(controlcard,"CASC").gt.0) then iz_sc=1 else if (index(controlcard,"SCONLY").gt.0) then iz_sc=2 else iz_sc=0 endif 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) pdbref=(index(controlcard,'PDBREF').gt.0) refstr = refstr .or. pdbref write (iout,*) "REFSTR",refstr," PDBREF",pdbref 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) print_fittest=(index(controlcard,'PRINT_FITTEST').gt.0) call readi(controlcard,'NCUT',ncut,0) if (ncut.gt.0) then call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) nclust=0 else call readi(controlcard,'NCLUST',nclust,5) endif 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) call readi(controlcard,'IOPT',iopt,2) lefree = 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_DIST',constr_dist,0) call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) c if (constr_homology) tole=dmax1(tole,1.5d0) write (iout,*) "with_homology_constr ",with_dihed_constr, & " CONSTR_HOMOLOGY",constr_homology read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD 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' include 'COMMON.SAXS' #ifdef MPL include 'COMMON.INFO' #endif character*4 sequence(maxres) character*800 weightcard,controlcard integer rescode double precision x(maxvar) double precision phihel,phibet,sigmahel,sigmabet,sumv, & secprob(3,maxres) integer itype_pdb(maxres) logical seq_comp integer i,j,kkk,i1,i2,it1,it2,tperm,ii,iperm C C Body C C Read weights of the subsequent energy terms. call card_concat(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) call reada(weightcard,'WDFAD',wdfa_dist,0.0d0) call reada(weightcard,'WDFAT',wdfa_tor,0.0d0) call reada(weightcard,'WDFAN',wdfa_nei,0.0d0) call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) 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) dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. enddo do i=1,max_cyst-1 do j=i+1,max_cyst 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 write (iout,*) "Parameters of the 'trisulfide' potential" write (iout,*) "ATRISS=", atriss write (iout,*) "BTRISS=", btriss write (iout,*) "CTRISS=", ctriss write (iout,*) "DTRISS=", dtriss 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)=scal14 weights(18)=wbond weights(19)=wsccor weights(28)=wdfa_dist weights(29)=wdfa_tor weights(30)=wdfa_nei weights(31)=wdfa_beta 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 #ifdef DFA write (iout,'(/a/)') "DFA pseudopotential parameters:" write (iout,'(a,f10.6,a)') & "WDFAD= ",wdfa_dist," (distance)", & "WDFAT= ",wdfa_tor," (backbone angles)", & "WDFAN= ",wdfa_nei," (neighbors)", & "WDFAB= ",wdfa_beta," (beta structure)" #endif call flush(iout) c 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 c print *,nres c 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 c print *,'Call Read_Bridge.' call read_bridge C this fragment reads diheadral constrains nnt=1 nct=nres c print *,'NNT=',NNT,' NCT=',NCT call seq2chains(nres,itype,nchain,chain_length,chain_border, & ireschain) write(iout,*) "nres",nres," nchain",nchain do i=1,nchain write(iout,*)"chain",i,chain_length(i),chain_border(1,i), & chain_border(2,i) enddo call chain_symmetry(nchain,nres,itype,chain_border, & chain_length,npermchain,tabpermchain) do i=1,nres write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain), & ii=1,npermchain) enddo write(iout,*) "residue permutations" do i=1,nres write(iout,*) i,(iperm(i,ii),ii=1,npermchain) enddo 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 nres0=nres #ifdef DFA if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and. & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then call init_dfa_vars print*, 'init_dfa_vars finished!' call read_dfa_info print*, 'read_dfa_info finished!' endif #endif C If the reference structure is not read set the superposition C boundaries nstart_sup=nnt nstart_seq=nnt nend_sup=nct nsup=nct-nnt+1 if (with_dihed_constr) then read (inp,*) ndih_constr if (ndih_constr.gt.0) then raw_psipred=.false. 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 else if (ndih_constr.lt.0) then raw_psipred=.true. call card_concat(controlcard) call reada(controlcard,"PHIHEL",phihel,50.0D0) call reada(controlcard,"PHIBET",phibet,180.0D0) call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0) call reada(controlcard,"SIGMABET",sigmabet,40.0d0) call reada(controlcard,"WDIHC",wdihc,0.591d0) write (iout,*) "Weight of the dihedral restraint term",wdihc read(inp,'(9x,3f7.3)') & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct) write (iout,*) "The secprob array" do i=nnt,nct write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3) enddo ndih_constr=0 do i=nnt+3,nct if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then ndih_constr=ndih_constr+1 idih_constr(ndih_constr)=i sumv=0.0d0 do j=1,3 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2) sumv=sumv+vpsipred(j,ndih_constr) enddo do j=1,3 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv enddo phibound(1,ndih_constr)=phihel*deg2rad phibound(2,ndih_constr)=phibet*deg2rad sdihed(1,ndih_constr)=sigmahel*deg2rad sdihed(2,ndih_constr)=sigmabet*deg2rad endif enddo write (iout,*) & 'There are',ndih_constr, & ' bimodal restraints on gamma angles.' do i=1,ndih_constr write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i, & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2, & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1, & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg, & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg, & (vpsipred(j,i),j=1,3) 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 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 if (constr_homology.gt.0) then c write (iout,*) "About to call read_constr_homology" c call flush(iout) call read_constr_homology c write (iout,*) "Exit read_constr_homology" c call flush(iout) if (indpdb.gt.0 .or. pdbref) then do i=1,2*nres do j=1,3 c(j,i)=crefjlee(j,i) cref(j,i)=crefjlee(j,i) enddo enddo endif #ifdef DEBUG write (iout,*) "Array C" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3), & (c(j,i+nres),j=1,3) enddo write (iout,*) "Array Cref" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3), & (cref(j,i+nres),j=1,3) enddo #endif #ifdef DEBUG call int_from_cart1(.false.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i) enddo do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) enddo enddo do i=2,nres-1 do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) enddo enddo #endif else homol_nset=0 endif write (iout,*) "calling read_saxs_consrtr",nsaxs if (nsaxs.gt.0) call read_saxs_constr 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 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) c print *,'ns=',ns c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine c numbers icys=0 do i=1,ns icys(iss(i))=i enddo 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') #ifndef NEWCORR call getenv('TORDPAR',tordname) open (itordp,file=tordname,status='old') #endif 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 c-------------------------------------------------------------------------- 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' include 'COMMON.INTERACT' 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 scal_bfac 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) call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0) if (restr_type.eq.10) & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0) if (restr_type.eq.12) & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0) 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 write (iout,*) "SCAL_PEAK",scal_peak write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ write (iout,*) "IFRAG" do i=1,nfrag_ write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) enddo write (iout,*) "IPAIR" do i=1,npair_ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) enddo if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) 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) 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)=wpair_(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)=wpair_(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)=wpair_(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.12) then read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), & fordepth_peak(nhpb_peak+1),npeak c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), c & fordepth_peak(nhpb_peak+1),npeak if (forcon_peak(nhpb_peak+1).le.0.0d0.or. & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle nhpb_peak=nhpb_peak+1 irestr_type_peak(nhpb_peak)=12 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i ipeak(2,npeak)=i write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) if (ibecarb_peak(nhpb_peak).eq.3) then jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres else if (ibecarb_peak(nhpb_peak).eq.2) then ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres else if (ibecarb_peak(nhpb_peak).eq.1) then ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres endif else 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) c if (ibecarb(nhpb).gt.0) then c ihpb(nhpb)=ihpb(nhpb)+nres c jhpb(nhpb)=jhpb(nhpb)+nres c endif if (ibecarb(nhpb).eq.3) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif else if (restr_type.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 c if (ibecarb(nhpb).gt.0) then c ihpb(nhpb)=ihpb(nhpb)+nres c jhpb(nhpb)=jhpb(nhpb)+nres c endif if (ibecarb(nhpb).eq.3) then jhpb(nhpb)=jhpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) 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 c if (ibecarb(nhpb).gt.0) then c ihpb(nhpb)=ihpb(nhpb)+nres c jhpb(nhpb)=jhpb(nhpb)+nres c endif if (ibecarb(nhpb).eq.3) then jhpb(nhpb)=jhpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) 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 if (restr_type.eq.4) then write (iout,*) "The BFAC array" do i=nnt,nct write (iout,'(i5,f10.5)') i,bfac(i) enddo do i=nnt,nct if (itype(i).eq.ntyp1) cycle do j=nnt,i-1 if (itype(j).eq.ntyp1) cycle if (itype(i).eq.10) then iiend=0 else iiend=1 endif if (itype(j).eq.10) then jjend=0 else jjend=1 endif kk=0 do ii=0,iiend do jj=0,jjend nhpb=nhpb+1 irestr_type(nhpb)=1 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2) irestr_type(nhpb)=1 ibecarb(nhpb)=kk if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb) ihpb(nhpb)=i+nres*ii jhpb(nhpb)=j+nres*jj dhpb(nhpb)=dist(i+nres*ii,j+nres*jj) 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) kk=kk+1 enddo enddo enddo enddo endif if (restr_type.eq.5) then restr_on_coord=.true. do i=nnt,nct if (itype(i).eq.ntyp1) cycle bfac(i)=(scal_bfac/bfac(i))**2 enddo endif 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' include 'COMMON.SAXS' 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