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