X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_Eshel%2Freadrtns.F;fp=source%2Funres%2Fsrc_Eshel%2Freadrtns.F;h=d321d987f564fcaca9ea89301fc003c76c152785;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_Eshel/readrtns.F b/source/unres/src_Eshel/readrtns.F new file mode 100644 index 0000000..d321d98 --- /dev/null +++ b/source/unres/src_Eshel/readrtns.F @@ -0,0 +1,1700 @@ + subroutine readrtns + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' + include 'COMMON.IOUNITS' + logical file_exist +C Read force-field parameters except weights + call parmread +C Read job setup parameters + call read_control +C Read control parameters for energy minimzation if required + if (minim) call read_minim +C Read molecule information, molecule geometry, energy-term weights, and +C restraints if requested + call molread +C Print restraint information +#ifdef MPI + if (.not. out1file .or. me.eq.king) then +#endif + if (nhpb.gt.nss) + &write (iout,'(a,i5,a)') "The following",nhpb-nss, + & " distance constraints have been imposed" + do i=nss+1,nhpb + write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i), + & ibecarb(i),dhpb(i),dhpb1(i),forcon(i) + enddo +#ifdef MPI + endif +#endif +c print *,"Processor",myrank," leaves READRTNS" + return + end +C------------------------------------------------------------------------------- + subroutine read_control +C +C Read contorl data +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MP + include 'mpif.h' + logical OKRandom, prng_restart + real*8 r1 +#endif + include 'COMMON.IOUNITS' + include 'COMMON.TIME1' + include 'COMMON.THREAD' + include 'COMMON.SBRIDGE' + include 'COMMON.CONTROL' + include 'COMMON.MCM' + include 'COMMON.MAP' + include 'COMMON.HEADER' +csa include 'COMMON.CSA' + include 'COMMON.CHAIN' + include 'COMMON.MUCA' + include 'COMMON.MD' + include 'COMMON.FFIELD' + include 'COMMON.SETUP' + COMMON /MACHSW/ KDIAG,ICORFL,IXDR + character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/ + character*80 ucase + character*320 controlcard + + nglob_csa=0 + eglob_csa=1d99 + nmin_csa=0 + read (INP,'(a)') titel + call card_concat(controlcard) +c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0 +c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file + call reada(controlcard,'SEED',seed,0.0D0) + call random_init(seed) +C Set up the time limit (caution! The time must be input in minutes!) + read_cart=index(controlcard,'READ_CART').gt.0 + catrace=index(controlcard,'CATRACE').gt.0 + call readi(controlcard,'CONSTR_DIST',constr_dist,0) + call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours + unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 + call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes + call reada(controlcard,'RMSDBC',rmsdbc,3.0D0) + call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0) + call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0) + call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0) + call reada(controlcard,'DRMS',drms,0.1D0) + if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc + write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 + write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max + write (iout,'(a,f10.1)')'DRMS = ',drms + write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm + write (iout,'(a,f10.1)') 'Time limit (min):',timlim + endif + call readi(controlcard,'NZ_START',nz_start,0) + call readi(controlcard,'NZ_END',nz_end,0) + call readi(controlcard,'IZ_SC',iz_sc,0) + timlim=60.0D0*timlim + safety = 60.0d0*safety + timem=timlim + modecalc=0 + call reada(controlcard,"T_BATH",t_bath,300.0d0) + minim=(index(controlcard,'MINIMIZE').gt.0) + dccart=(index(controlcard,'CART').gt.0) + overlapsc=(index(controlcard,'OVERLAP').gt.0) + overlapsc=.not.overlapsc + searchsc=(index(controlcard,'SEARCHSC').gt.0) + sideadd=(index(controlcard,'SIDEADD').gt.0) .or. catrace + energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) + outpdb=(index(controlcard,'PDBOUT').gt.0) + outmol2=(index(controlcard,'MOL2OUT').gt.0) + pdbref=(index(controlcard,'PDBREF').gt.0) + refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0) + indpdb=index(controlcard,'PDBSTART') + extconf=(index(controlcard,'EXTCONF').gt.0) + call readi(controlcard,'IPRINT',iprint,0) + call readi(controlcard,'MAXGEN',maxgen,10000) + call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) + call readi(controlcard,"RESCALE_MODE",rescale_mode,2) + if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) + & write (iout,*) "RESCALE_MODE",rescale_mode + if (index(controlcard,'REGULAR').gt.0.0D0) then + call reada(controlcard,'WEIDIS',weidis,0.1D0) + modecalc=1 + refstr=.true. + regular=.true. + endif + if (index(controlcard,'CHECKGRAD').gt.0) then + modecalc=5 + if (index(controlcard,'CART').gt.0) then + icheckgrad=1 + elseif (index(controlcard,'CARINT').gt.0) then + icheckgrad=2 + else + icheckgrad=3 + endif + elseif (index(controlcard,'THREAD').gt.0) then + modecalc=2 + call readi(controlcard,'THREAD',nthread,0) + if (nthread.gt.0) then + call reada(controlcard,'WEIDIS',weidis,0.1D0) + else + if (fg_rank.eq.0) + & write (iout,'(a)')'A number has to follow the THREAD keyword.' + stop 'Error termination in Read_Control.' + endif + else if (index(controlcard,'MCMA').gt.0) then + modecalc=3 + else if (index(controlcard,'MCEE').gt.0) then + modecalc=6 + else if (index(controlcard,'MULTCONF').gt.0) then + modecalc=4 + else if (index(controlcard,'MAP').gt.0) then + modecalc=7 + call readi(controlcard,'MAP',nmap,0) + else if (index(controlcard,'CSA').gt.0) then + write(*,*) "CSA not supported in this version" + stop + else if (index(controlcard,'SOFTREG').gt.0) then + modecalc=11 + else if (index(controlcard,'CHECK_BOND').gt.0) then + modecalc=-1 + else if (index(controlcard,'TEST').gt.0) then + modecalc=-2 + endif + + lmuca=index(controlcard,'MUCA').gt.0 + call readi(controlcard,'MUCADYN',mucadyn,0) + call readi(controlcard,'MUCASMOOTH',muca_smooth,0) + if (lmuca .and. (me.eq.king .or. .not.out1file )) + & then + write (iout,*) 'MUCADYN=',mucadyn + write (iout,*) 'MUCASMOOTH=',muca_smooth + endif + + iscode=index(controlcard,'ONE_LETTER') + indphi=index(controlcard,'PHI') + indback=index(controlcard,'BACK') + iranconf=index(controlcard,'RAND_CONF') + i2ndstr=index(controlcard,'USE_SEC_PRED') + gradout=index(controlcard,'GRADOUT').gt.0 + gnorm_check=index(controlcard,'GNORM_CHECK').gt.0 + + if(me.eq.king.or..not.out1file) + & write (iout,'(2a)') diagmeth(kdiag), + & ' routine used to diagonalize matrices.' + return + end +c-------------------------------------------------------------------------- + subroutine molread +C +C Read molecular data. +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' + integer error_msg +#endif + 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.DBASE' + include 'COMMON.THREAD' + include 'COMMON.CONTACTS' + include 'COMMON.TORCNSTR' + include 'COMMON.TIME1' + include 'COMMON.BOUNDS' + include 'COMMON.MD' + include 'COMMON.REMD' + include 'COMMON.SETUP' + character*4 sequence(maxres) + integer rescode + double precision x(maxvar) + character*320 weightcard + character*80 weightcard_t,ucase + dimension itype_pdb(maxres) + common /pizda/ itype_pdb + logical seq_comp,fail + double precision energia(0:n_ene) + integer ilen + external ilen +C +C Body +C +C Read weights of the subsequent energy terms. + if(hremd.gt.0) then + + k=0 + do il=1,hremd + do i=1,nrep + do j=1,remd_m(i) + i2set(k)=il + k=k+1 + enddo + enddo + enddo + + if(me.eq.king.or..not.out1file) then + write (iout,*) 'Reading ',hremd,' sets of weights for HREMD' + write (iout,*) 'Current weights for processor ', + & me,' set ',i2set(me) + endif + + do i=1,hremd + call card_concat(weightcard) + call reada(weightcard,'WLONG',wlong,1.0D0) + call reada(weightcard,'WSC',wsc,wlong) + call reada(weightcard,'WSCP',wscp,wlong) + 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,'WSCCOR',wsccor,1.0D0) + call reada(weightcard,'WSTRAIN',wstrain,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,'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) + call reada(weightcard,'TEMP0',temp0,300.0d0) + if (index(weightcard,'SOFT').gt.0) ipot=6 +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 + + hweights(i,1)=wsc + hweights(i,2)=wscp + hweights(i,3)=welec + hweights(i,4)=wcorr + hweights(i,5)=wcorr5 + hweights(i,6)=wcorr6 + hweights(i,7)=wel_loc + hweights(i,8)=wturn3 + hweights(i,9)=wturn4 + hweights(i,10)=wturn6 + hweights(i,11)=wang + hweights(i,12)=wscloc + hweights(i,13)=wtor + hweights(i,14)=wtor_d + hweights(i,15)=wstrain + hweights(i,16)=wvdwpp + hweights(i,17)=wbond + hweights(i,18)=scal14 + hweights(i,21)=wsccor + + enddo + + do i=1,n_ene + weights(i)=hweights(i2set(me),i) + enddo + wsc =weights(1) + wscp =weights(2) + welec =weights(3) + wcorr =weights(4) + wcorr5 =weights(5) + wcorr6 =weights(6) + wel_loc=weights(7) + wturn3 =weights(8) + wturn4 =weights(9) + wturn6 =weights(10) + wang =weights(11) + wscloc =weights(12) + wtor =weights(13) + wtor_d =weights(14) + wstrain=weights(15) + wvdwpp =weights(16) + wbond =weights(17) + scal14 =weights(18) + wsccor =weights(21) + + + else + call card_concat(weightcard) + call reada(weightcard,'WLONG',wlong,1.0D0) + call reada(weightcard,'WSC',wsc,wlong) + call reada(weightcard,'WSCP',wscp,wlong) + 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,'WSCCOR',wsccor,1.0D0) + call reada(weightcard,'WSTRAIN',wstrain,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,'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) + call reada(weightcard,'TEMP0',temp0,300.0d0) + if (index(weightcard,'SOFT').gt.0) ipot=6 +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(21)=wsccor + endif + + if(me.eq.king.or..not.out1file) + & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, + & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, + & wturn4,wturn6 + 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)'/ + & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ + & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ + & 'WTURN4= ',f10.6,' (turns, 4th order)'/ + & 'WTURN6= ',f10.6,' (turns, 6th order)') + if(me.eq.king.or..not.out1file)then + 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 + endif + 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 rescale_weights(t_bath) + if(me.eq.king.or..not.out1file) + & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, + & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, + & wturn4,wturn6 + 22 format (/'Energy-term weights (scaled):'// + & '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)'/ + & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ + & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ + & 'WTURN4= ',f10.6,' (turns, 4th order)'/ + & 'WTURN6= ',f10.6,' (turns, 6th order)') + if(me.eq.king.or..not.out1file) + & write (iout,*) "Reference temperature for weights calculation:", + & temp0 + 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) + if(me.eq.king.or..not.out1file) then + write (iout,*) "Parameters of the SS-bond potential:" + write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, + & " AKCT",akct + write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss + write (iout,*) "EBR",ebr + print *,'indpdb=',indpdb,' pdbref=',pdbref + endif + if (indpdb.gt.0 .or. pdbref) then + read(inp,'(a)') pdbname(1) + call pdbinput(1) + call contact(.false.,ncont_ref,icont_ref,co) + endif + if (indpdb.eq.0) then +C Read sequence if not taken from the pdb file. + read (inp,*) nres +c print *,'nres=',nres + 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 Assign initial virtual bond lengths + do i=2,nres + vbld(i)=vbl + vbld_inv(i)=vblinv + enddo + do i=2,nres-1 + vbld(i+nres)=dsc(itype(i)) + vbld_inv(i+nres)=dsc_inv(itype(i)) +c write (iout,*) "i",i," itype",itype(i), +c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres) + enddo + endif +c print *,nres +c print '(20i4)',(itype(i),i=1,nres) + do i=1,nres +#ifdef PROCOR + if (itype(i).eq.21 .or. itype(i+1).eq.21) then +#else + if (itype(i).eq.21) then +#endif + itel(i)=0 +#ifdef PROCOR + else if (itype(i+1).ne.20) then +#else + else if (itype(i).ne.20) then +#endif + itel(i)=1 + else + itel(i)=2 + endif + enddo + if(me.eq.king.or..not.out1file)then + write (iout,*) "ITEL" + do i=1,nres-1 + write (iout,*) i,itype(i),itel(i) + enddo + print *,'Call Read_Bridge.' + endif + call read_bridge +C 8/13/98 Set limits to generating the dihedral angles + do i=1,nres + phibound(1,i)=-pi + phibound(2,i)=pi + enddo + read (inp,*) ndih_constr + if (ndih_constr.gt.0) then + read (inp,*) ftors + read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) + if(me.eq.king.or..not.out1file)then + write (iout,*) + & 'There are',ndih_constr,' constraints on phi angles.' + do i=1,ndih_constr + write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) + enddo + endif + do i=1,ndih_constr + phi0(i)=deg2rad*phi0(i) + drange(i)=deg2rad*drange(i) + enddo + if(me.eq.king.or..not.out1file) + & write (iout,*) 'FTORS',ftors + do i=1,ndih_constr + ii = idih_constr(i) + phibound(1,ii) = phi0(i)-drange(i) + phibound(2,ii) = phi0(i)+drange(i) + enddo + endif + nnt=1 +#ifdef MPI + if (me.eq.king) then +#endif + write (iout,'(a)') 'Boundaries in phi angle sampling:' + do i=1,nres + write (iout,'(a3,i5,2f10.1)') + & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg + enddo +#ifdef MP + endif +#endif + nct=nres +cd print *,'NNT=',NNT,' NCT=',NCT + if (itype(1).eq.21) nnt=2 + if (itype(nres).eq.21) nct=nct-1 + if (pdbref) then + if(me.eq.king.or..not.out1file) + & write (iout,'(a,i3)') 'nsup=',nsup + nstart_seq=nnt + if (nsup.le.(nct-nnt+1)) then + do i=0,nct-nnt+1-nsup + if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then + nstart_seq=nnt+i + goto 111 + endif + enddo + write (iout,'(a)') + & 'Error - sequences to be superposed do not match.' + stop + else + do i=0,nsup-(nct-nnt+1) + if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) + & then + nstart_sup=nstart_sup+i + nsup=nct-nnt+1 + goto 111 + endif + enddo + write (iout,'(a)') + & 'Error - sequences to be superposed do not match.' + endif + 111 continue + if (nsup.eq.0) nsup=nct-nnt + if (nstart_sup.eq.0) nstart_sup=nnt + if (nstart_seq.eq.0) nstart_seq=nnt + if(me.eq.king.or..not.out1file) + & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup, + & ' nstart_seq=',nstart_seq + endif +c--- Zscore rms ------- + if (nz_start.eq.0) nz_start=nnt + if (nz_end.eq.0 .and. nsup.gt.0) then + nz_end=nnt+nsup-1 + else if (nz_end.eq.0) then + nz_end=nct + endif + if(me.eq.king.or..not.out1file)then + write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end + write (iout,*) 'IZ_SC=',iz_sc + endif +c---------------------- + call init_int_table + if (refstr) then + if (.not.pdbref) then + call read_angles(inp,*38) + goto 39 + 38 write (iout,'(a)') 'Error reading reference structure.' +#ifdef MPI + call MPI_Finalize(MPI_COMM_WORLD,IERROR) + stop 'Error reading reference structure' +#endif + 39 call chainbuild + call setup_var +czscore call geom_to_var(nvar,coord_exp_zs(1,1)) + nstart_sup=nnt + nstart_seq=nnt + nsup=nct-nnt+1 + do i=1,2*nres + do j=1,3 + cref(j,i)=c(j,i) + enddo + enddo + call contact(.true.,ncont_ref,icont_ref,co) + endif + if(me.eq.king.or..not.out1file) + & write (iout,*) 'Contact order:',co + if (pdbref) then + if(me.eq.king.or..not.out1file) + & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup + do i=1,ncont_ref + do j=1,2 + icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup + enddo + if(me.eq.king.or..not.out1file) + & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ', + & icont_ref(1,i),' ', + & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i) + enddo + endif + endif +c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup + if (constr_dist.gt.0) then + call read_dist_constr + endif + if (nhpb.gt.0) call hpb_partition +c write (iout,*) "After read_dist_constr nhpb",nhpb +c call flush(iout) + if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 + & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. + & modecalc.ne.10) then +C If input structure hasn't been supplied from the PDB file read or generate +C initial geometry. + if (iranconf.eq.0 .and. .not. extconf) then + if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) + & write (iout,'(a)') 'Initial geometry will be read in.' + if (read_cart) then + read(inp,'(8f10.5)',end=36,err=36) + & ((c(l,k),l=1,3),k=1,nres), + & ((c(l,k+nres),l=1,3),k=nnt,nct) + call int_from_cart1(.false.) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1) + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres) + enddo + endif + enddo + return + else + call read_angles(inp,*36) + endif + goto 37 + 36 write (iout,'(a)') 'Error reading angle file.' +#ifdef MPI + call mpi_finalize( MPI_COMM_WORLD,IERR ) +#endif + stop 'Error reading angle file.' + 37 continue + else if (extconf) then + if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) + & write (iout,'(a)') 'Extended chain initial geometry.' + do i=3,nres + theta(i)=90d0*deg2rad + enddo + do i=4,nres + phi(i)=180d0*deg2rad + enddo + do i=2,nres-1 + alph(i)=110d0*deg2rad + enddo + do i=2,nres-1 + omeg(i)=-120d0*deg2rad + enddo + else + if(me.eq.king.or..not.out1file) + & write (iout,'(a)') 'Random-generated initial geometry.' + + +#ifdef MPI + if (me.eq.king .or. fg_rank.eq.0 .and. ( + & modecalc.eq.12 .or. modecalc.eq.14) ) then +#endif + do itrial=1,100 + itmp=1 + call gen_rand_conf(itmp,*30) + goto 40 + 30 write (iout,*) 'Failed to generate random conformation', + & ', itrial=',itrial + write (*,*) 'Processor:',me, + & ' Failed to generate random conformation', + & ' itrial=',itrial + call intout + +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif + enddo + write (iout,'(a,i3,a)') 'Processor:',me, + & ' error in generating random conformation.' + write (*,'(a,i3,a)') 'Processor:',me, + & ' error in generating random conformation.' + call flush(iout) +#ifdef MPI + call MPI_Abort(mpi_comm_world,error_msg,ierrcode) + 40 continue + endif +#else + 40 continue +#endif + endif + elseif (modecalc.eq.4) then + read (inp,'(a)') intinname + open (intin,file=intinname,status='old',err=333) + if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) + & write (iout,'(a)') 'intinname',intinname + write (*,'(a)') 'Processor',myrank,' intinname',intinname + goto 334 + 333 write (iout,'(2a)') 'Error opening angle file ',intinname +#ifdef MPI + call MPI_Finalize(MPI_COMM_WORLD,IERR) +#endif + stop 'Error opening angle file.' + 334 continue + + endif +C Generate distance constraints, if the PDB structure is to be regularized. + if (nthread.gt.0) then + call read_threadbase + endif + call setup_var + if (me.eq.king .or. .not. out1file) + & call intout + if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then + write (iout,'(/a,i3,a)') + & 'The chain contains',ns,' disulfide-bridging cysteines.' + write (iout,'(20i4)') (iss(i),i=1,ns) + 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) + if (me.eq.king.or..not.out1file) + & 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 + if (i2ndstr.gt.0) call secstrp2dihc +c call geom_to_var(nvar,x) +c call etotal(energia(0)) +c call enerprint(energia(0)) +c call briefout(0,etot) +c stop +cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT +cd write (iout,'(a)') 'Variable list:' +cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar) + if (indpdb.gt.0) then + do i=1,10000 + read(1,'(a)',end=555,err=555) pdbname(i) + enddo + 555 npdbfile=i-1 + write (iout,*) "PDB files" + do i=1,npdbfile + write (iout,'(i6,a)') i,pdbname(i)(:ilen(pdbname(i))) + enddo + endif +#ifdef MPI + if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) + & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') + & 'Processor',myrank,': end reading molecular data.' +#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 real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + 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.DBASE' + include 'COMMON.THREAD' + include 'COMMON.TIME1' + include 'COMMON.SETUP' +C Read bridging residues. + read (inp,*) ns,(iss(i),i=1,ns) + print *,'ns=',ns + if(me.eq.king.or..not.out1file) + & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns) +C Check whether the specified bridging residues are cystines. + do i=1,ns + if (itype(iss(i)).ne.1) then + if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') + & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & ' can form a disulfide bridge?!!!' + write (*,'(2a,i3,a)') + & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & ' can form a disulfide bridge?!!!' +#ifdef MPI + call MPI_Finalize(MPI_COMM_WORLD,ierror) + stop +#endif + endif + enddo +C Read preformed bridges. + if (ns.gt.0) then + read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) + write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(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 MPI + call MPI_Finalize(MPI_COMM_WORLD,ierror) + 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 + dhpb(i)=dbr + 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_x(kanal,*) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' +c Read coordinates from input +c + read(kanal,'(8f10.5)',end=10,err=10) + & ((c(l,k),l=1,3),k=1,nres), + & ((c(l,k+nres),l=1,3),k=nnt,nct) + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + call int_from_cart1(.false.) + 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=nnt,nct + if (itype(i).ne.10) then + 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 + endif + enddo + + return + 10 return1 + end +c---------------------------------------------------------------------------- + subroutine read_threadbase + implicit real*8 (a-h,o-z) + 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.DBASE' + include 'COMMON.THREAD' + include 'COMMON.TIME1' +C Read pattern database for threading. + read (icbase,*) nseq + do i=1,nseq + read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i), + & nres_base(2,i),nres_base(3,i) + read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1, + & nres_base(1,i)) +c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i), +c & nres_base(2,i),nres_base(3,i) +c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1, +c & nres_base(1,i)) + enddo + close (icbase) + if (weidis.eq.0.0D0) weidis=0.1D0 + do i=nnt,nct + do j=i+2,nct + nhpb=nhpb+1 + ihpb(nhpb)=i + jhpb(nhpb)=j + forcon(nhpb)=weidis + enddo + enddo + read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl) + write (iout,'(a,i5)') 'nexcl: ',nexcl + write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl) + return + end +c------------------------------------------------------------------------------ + subroutine setup_var + implicit real*8 (a-h,o-z) + 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.DBASE' + include 'COMMON.THREAD' + include 'COMMON.TIME1' +C Set up variable list. + ntheta=nres-2 + nphi=nres-3 + nvar=ntheta+nphi + nside=0 + do i=2,nres-1 + if (itype(i).ne.10) then + nside=nside+1 + ialph(i,1)=nvar+nside + ialph(nside,2)=i + endif + enddo + if (indphi.gt.0) then + nvar=nphi + else if (indback.gt.0) then + nvar=nphi+ntheta + else + nvar=nvar+2*nside + endif +cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1) + return + end +c---------------------------------------------------------------------------- + subroutine gen_dist_constr +C Generate CA distance constraints. + implicit real*8 (a-h,o-z) + 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.DBASE' + include 'COMMON.THREAD' + include 'COMMON.TIME1' + dimension itype_pdb(maxres) + common /pizda/ itype_pdb + character*2 iden +cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct +cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct, +cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq, +cd & ' nsup',nsup + do i=nstart_sup,nstart_sup+nsup-1 +cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)), +cd & ' seq_pdb', restyp(itype_pdb(i)) + do j=i+2,nstart_sup+nsup-1 + nhpb=nhpb+1 + ihpb(nhpb)=i+nstart_seq-nstart_sup + jhpb(nhpb)=j+nstart_seq-nstart_sup + forcon(nhpb)=weidis + dhpb(nhpb)=dist(i,j) + enddo + enddo +cd write (iout,'(a)') 'Distance constraints:' +cd do i=nss+1,nhpb +cd ii=ihpb(i) +cd jj=jhpb(i) +cd iden='CA' +cd if (ii.gt.nres) then +cd iden='SC' +cd ii=ii-nres +cd jj=jj-nres +cd endif +cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)') +cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj, +cd & dhpb(i),forcon(i) +cd enddo + return + end +c---------------------------------------------------------------------------- + subroutine map_read + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.MAP' + include 'COMMON.IOUNITS' + character*3 angid(4) /'THE','PHI','ALP','OME'/ + character*80 mapcard,ucase + do imap=1,nmap + read (inp,'(a)') mapcard + mapcard=ucase(mapcard) + if (index(mapcard,'PHI').gt.0) then + kang(imap)=1 + else if (index(mapcard,'THE').gt.0) then + kang(imap)=2 + else if (index(mapcard,'ALP').gt.0) then + kang(imap)=3 + else if (index(mapcard,'OME').gt.0) then + kang(imap)=4 + else + write(iout,'(a)')'Error - illegal variable spec in MAP card.' + stop 'Error - illegal variable spec in MAP card.' + endif + call readi (mapcard,'RES1',res1(imap),0) + call readi (mapcard,'RES2',res2(imap),0) + if (res1(imap).eq.0) then + res1(imap)=res2(imap) + else if (res2(imap).eq.0) then + res2(imap)=res1(imap) + endif + if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then + write (iout,'(a)') + & 'Error - illegal definition of variable group in MAP.' + stop 'Error - illegal definition of variable group in MAP.' + endif + call reada(mapcard,'FROM',ang_from(imap),0.0D0) + call reada(mapcard,'TO',ang_to(imap),0.0D0) + call readi(mapcard,'NSTEP',nstep(imap),0) + if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then + write (iout,'(a)') + & 'Illegal boundary and/or step size specification in MAP.' + stop 'Illegal boundary and/or step size specification in MAP.' + endif + enddo ! imap + return + end +c---------------------------------------------------------------------------- + subroutine read_minim + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.MINIM' + include 'COMMON.IOUNITS' + character*80 ucase + character*320 minimcard + call card_concat(minimcard) + call readi(minimcard,'MAXMIN',maxmin,2000) + call readi(minimcard,'MAXFUN',maxfun,5000) + call readi(minimcard,'MINMIN',minmin,maxmin) + call readi(minimcard,'MINFUN',minfun,maxmin) + call reada(minimcard,'TOLF',tolf,1.0D-2) + call reada(minimcard,'RTOLF',rtolf,1.0D-4) + print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1) + print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1) + print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1) + write (iout,'(/80(1h*)/20x,a/80(1h*))') + & 'Options in energy minimization:' + write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') + & 'MaxMin:',MaxMin,' MaxFun:',MaxFun, + & 'MinMin:',MinMin,' MinFun:',MinFun, + & ' TolF:',TolF,' RTolF:',RTolF + return + end +c---------------------------------------------------------------------------- + subroutine read_angles(kanal,*) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' +c Read angles from input +c + 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 +c 9/7/01 avoid 180 deg valence angle + if (theta(i).gt.179.99d0) theta(i)=179.99d0 +c + 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:),*,err=10,end=10) wartosc + return + 10 wartosc=default + 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:),*,err=10,end=10) wartosc + return + 10 wartosc=default + 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 multreada(rekord,lancuch,tablica,dim,default) + implicit none + integer dim,i + double precision 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) + implicit real*8 (a-h,o-z) + 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 read_fragments + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.CONTROL' + read(inp,*) nset,nfrag,npair,nfrag_back + if(me.eq.king.or..not.out1file) + & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair, + & " nfrag_back",nfrag_back + do iset=1,nset + read(inp,*) mset(iset) + do i=1,nfrag + read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset), + & qinfrag(i,iset) + if(me.eq.king.or..not.out1file) + & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset), + & ifrag(2,i,iset), qinfrag(i,iset) + enddo + do i=1,npair + read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset), + & qinpair(i,iset) + if(me.eq.king.or..not.out1file) + & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset), + & ipair(2,i,iset), qinpair(i,iset) + enddo + do i=1,nfrag_back + read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset), + & wfrag_back(3,i,iset), + & ifrag_back(1,i,iset),ifrag_back(2,i,iset) + if(me.eq.king.or..not.out1file) + & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset), + & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset) + enddo + enddo + return + end +c------------------------------------------------------------------------------- + subroutine read_dist_constr + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + 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 +c write (iout,*) "Calling read_dist_constr" +c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +c call flush(iout) + call card_concat(controlcard) + 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 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) +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 (.not.refstr .and. nfrag.gt.0) then + write (iout,*) + & "ERROR: no reference structure to compute distance restraints" + write (iout,*) + & "Restraints must be specified explicitly (NDIST=number)" + stop + endif + if (nfrag.lt.2 .and. npair.gt.0) then + write (iout,*) "ERROR: Less than 2 fragments specified", + & " but distance restraints between pairs requested" + stop + endif + call flush(iout) + 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) + call flush(iout) + if (wfrag_(i).gt.0.0d0) then + do j=ifrag_(1,i),ifrag_(2,i)-1 + do k=j+1,ifrag_(2,i) + write (iout,*) "j",j," k",k + ddjk=dist(j,k) + if (constr_dist.eq.1) then + nhpb=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 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + endif + else + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) + endif +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#else + write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#endif + enddo + enddo + endif + enddo + do i=1,npair_ + if (wpair_(i).gt.0.0d0) then + 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) + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + forcon(nhpb)=wpair_(i) + dhpb(nhpb)=dist(j,k) +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#else + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#endif + enddo + enddo + endif + enddo + do i=1,ndist_ + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) + if (forcon(nhpb+1).gt.0.0d0) then + nhpb=nhpb+1 + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif + enddo +#ifdef MPI + if (.not.out1file .or. me.eq.king) then +#endif + do i=1,nhpb + write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", + & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) + enddo + call flush(iout) +#ifdef MPI + endif +#endif + return + end +c------------------------------------------------------------------------------- +#ifdef WINIFL + subroutine flush(iu) + return + end +#endif +#ifdef AIX + subroutine flush(iu) + call flush_(iu) + return + end +#endif +c------------------------------------------------------------------------------ + subroutine copy_to_tmp(source) + include "DIMENSIONS" + include "COMMON.IOUNITS" + character*(*) source + character* 256 tmpfile + integer ilen + external ilen + logical ex + tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source)) + inquire(file=tmpfile,exist=ex) + if (ex) then + write (*,*) "Copying ",tmpfile(:ilen(tmpfile)), + & " to temporary directory..." + write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir + call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir) + endif + return + end +c------------------------------------------------------------------------------ + subroutine move_from_tmp(source) + include "DIMENSIONS" + include "COMMON.IOUNITS" + character*(*) source + integer ilen + external ilen + write (*,*) "Moving ",source(:ilen(source)), + & " from temporary directory to working directory" + write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir + call system("/bin/mv "//source(:ilen(source))//" "//curdir) + return + end +c------------------------------------------------------------------------------ + subroutine random_init(seed) +C +C Initialize random number generator +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef AMD64 + integer*8 iseedi8 +#endif +#ifdef MPI + include 'mpif.h' + logical OKRandom, prng_restart + real*8 r1 + integer iseed_array(4) +#endif + include 'COMMON.IOUNITS' + include 'COMMON.TIME1' + include 'COMMON.THREAD' + include 'COMMON.SBRIDGE' + include 'COMMON.CONTROL' + include 'COMMON.MCM' + include 'COMMON.MAP' + include 'COMMON.HEADER' +csa include 'COMMON.CSA' + include 'COMMON.CHAIN' + include 'COMMON.MUCA' + include 'COMMON.MD' + include 'COMMON.FFIELD' + include 'COMMON.SETUP' + iseed=-dint(dabs(seed)) + if (iseed.eq.0) then + write (iout,'(/80(1h*)/20x,a/80(1h*))') + & 'Random seed undefined. The program will stop.' + write (*,'(/80(1h*)/20x,a/80(1h*))') + & 'Random seed undefined. The program will stop.' +#ifdef MPI + call mpi_finalize(mpi_comm_world,ierr) +#endif + stop 'Bad random seed.' + endif +#ifdef MPI + if (fg_rank.eq.0) then + seed=seed*(me+1)+1 +#ifdef AMD64 + iseedi8=dint(seed) + if(me.eq.king .or. .not. out1file) + & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8 + write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8 + OKRandom = prng_restart(me,iseedi8) +#else + do i=1,4 + tmp=65536.0d0**(4-i) + iseed_array(i) = dint(seed/tmp) + seed=seed-iseed_array(i)*tmp + enddo + if(me.eq.king .or. .not. out1file) + & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ', + & (iseed_array(i),i=1,4) + write (*,*) 'MPI: node= ',me, ' iseed(4)= ', + & (iseed_array(i),i=1,4) + OKRandom = prng_restart(me,iseed_array) +#endif + if (OKRandom) then + r1=ran_number(0.0D0,1.0D0) + if(me.eq.king .or. .not. out1file) + & write (iout,*) 'ran_num',r1 + if (r1.lt.0.0d0) OKRandom=.false. + endif + if (.not.OKRandom) then + write (iout,*) 'PRNG IS NOT WORKING!!!' + print *,'PRNG IS NOT WORKING!!!' + if (me.eq.0) then + call flush(iout) + call mpi_abort(mpi_comm_world,error_msg,ierr) + stop + else + write (iout,*) 'too many processors for parallel prng' + write (*,*) 'too many processors for parallel prng' + call flush(iout) + stop + endif + endif + endif +#else + call vrndst(iseed) + write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0) +#endif + return + end +c-------------------------------------------------------------------------- + subroutine pdbinput(iconf) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' + integer error_msg +#endif + 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.DBASE' + include 'COMMON.THREAD' + include 'COMMON.CONTACTS' + include 'COMMON.TORCNSTR' + include 'COMMON.TIME1' + include 'COMMON.BOUNDS' + include 'COMMON.MD' + include 'COMMON.REMD' + include 'COMMON.SETUP' + character*4 sequence(maxres) + integer rescode + double precision x(maxvar) + dimension itype_pdb(maxres) + common /pizda/ itype_pdb + logical seq_comp,fail + double precision energia(0:n_ene) + integer ilen + external ilen + if(me.eq.king.or..not.out1file) + & write (iout,'(2a)') 'PDB data will be read from file ', + & pdbname(iconf)(:ilen(pdbname(iconf))) + open(ipdbin,file=pdbname(iconf),status='old',err=33) + goto 34 + 33 write (iout,'(a)') 'Error opening PDB file.' + stop + 34 continue +c print *,'Begin reading pdb data' + call readpdb +c print *,'Finished reading pdb data' + if(me.eq.king.or..not.out1file) + & write (iout,'(a,i3,a,i3)')'nsup=',nsup, + & ' nstart_sup=',nstart_sup + do i=1,nres + itype_pdb(i)=itype(i) + enddo + close (ipdbin) + nnt=nstart_sup + nct=nstart_sup+nsup-1 + + if (sideadd) then +C Following 2 lines for diagnostics; comment out if not needed + write (iout,*) "Before sideadd" + call intout + if (me.eq.king.or..not.out1file) + & write(iout,*)'Adding sidechains' + maxsi=1000 + do i=2,nres-1 + iti=itype(i) + if (iti.ne.10) then + nsi=0 + fail=.true. + do while (fail.and.nsi.le.maxsi) + call gen_side(iti,theta(i+1),alph(i),omeg(i),fail) + nsi=nsi+1 + enddo + if(fail) write(iout,*)'Adding sidechain failed for res ', + & i,' after ',nsi,' trials' + endif + enddo +C 10/03/12 Adam: Recalculate coordinates with new side chain positions + call chainbuild +C Following 2 lines for diagnostics; comment out if not needed + write (iout,*) "After sideadd" + call intout + endif + return + end