--- /dev/null
+ 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
projects / unres.git / blobdiff