include 'COMMON.CONTROL'
include 'COMMON.SBRIDGE'
include 'COMMON.IOUNITS'
+ include 'COMMON.SPLITELE'
logical file_exist
-C Read force-field parameters except weights
- call parmread
C Read job setup parameters
call read_control
+C Read force-field parameters except weights
+ call parmread
C Read control parameters for energy minimzation if required
if (minim) call read_minim
C Read MCM control parameters if required
include 'COMMON.FFIELD'
include 'COMMON.INTERACT'
include 'COMMON.SETUP'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SHIELD'
+ include 'COMMON.GEO'
COMMON /MACHSW/ KDIAG,ICORFL,IXDR
character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
character*80 ucase
C Set up the time limit (caution! The time must be input in minutes!)
read_cart=index(controlcard,'READ_CART').gt.0
call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+C this variable with_theta_constr is the variable which allow to read and execute the
+C constrains on theta angles WITH_THETA_CONSTR is the keyword
+ with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+ write (iout,*) "with_theta_constr ",with_theta_constr
call readi(controlcard,'SYM',symetr,1)
- call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
+ call reada(controlcard,'TIMLIM',timlim,2800.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)
searchsc=.not.searchsc
sideadd=(index(controlcard,'SIDEADD').gt.0)
energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
+ mremd_dec=(index(controlcard,'MREMD_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')
+ call readi(controlcard,'PDBSTART',indpdb,0)
+ if (index(controlcard,'PDBSTART').gt.0 .and. indpdb.eq.0) indpdb=2
+ write (iout,*) "indpdb",indpdb
extconf=(index(controlcard,'EXTCONF').gt.0)
+ AFMlog=(index(controlcard,'AFM'))
+ selfguide=(index(controlcard,'SELFGUIDE'))
+c print *,'AFMlog',AFMlog,selfguide,"KUPA"
+ call readi(controlcard,'TUBEMOD',tubelog,0)
+c write (iout,*) TUBElog,"TUBEMODE"
call readi(controlcard,'IPRINT',iprint,0)
+C SHIELD keyword sets if the shielding effect of side-chains is used
+C 0 denots no shielding is used all peptide are equally despite the
+C solvent accesible area
+C 1 the newly introduced function
+C 2 reseved for further possible developement
+ call readi(controlcard,'SHIELD',shield_mode,0)
+C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "shield_mode",shield_mode
+C endif
+ call readi(controlcard,'TORMODE',tor_mode,0)
+C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "torsional and valence angle mode",tor_mode
call readi(controlcard,'MAXGEN',maxgen,10000)
call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
call readi(controlcard,"KDIAG",kdiag,0)
endif
if (index(controlcard,'CHECKGRAD').gt.0) then
modecalc=5
- if (index(controlcard,'CART').gt.0) then
+ if (index(controlcard,' CART').gt.0) then
icheckgrad=1
elseif (index(controlcard,'CARINT').gt.0) then
icheckgrad=2
else
icheckgrad=3
endif
+ call reada(controlcard,'DELTA',aincr,1.0d-7)
+c write (iout,*) "icheckgrad",icheckgrad
elseif (index(controlcard,'THREAD').gt.0) then
modecalc=2
call readi(controlcard,'THREAD',nthread,0)
i2ndstr=index(controlcard,'USE_SEC_PRED')
gradout=index(controlcard,'GRADOUT').gt.0
gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
+C DISTCHAINMAX become obsolete for periodic boundry condition
call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
+C Reading the dimensions of box in x,y,z coordinates
+ call reada(controlcard,'BOXX',boxxsize,100.0d0)
+ call reada(controlcard,'BOXY',boxysize,100.0d0)
+ call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+c Cutoff range for interactions
+ call reada(controlcard,"R_CUT",r_cut,15.0d0)
+ call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+ call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+ call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+ if (lipthick.gt.0.0d0) then
+ bordliptop=(boxzsize+lipthick)/2.0
+ bordlipbot=bordliptop-lipthick
+C endif
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
+ & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+ buflipbot=bordlipbot+lipbufthick
+ bufliptop=bordliptop-lipbufthick
+ if ((lipbufthick*2.0d0).gt.lipthick)
+ &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+ endif
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+ write (iout,*) "SHIELD MODE",shield_mode
+ if (TUBElog.gt.0) then
+ call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+ call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
+ call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+ call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
+ call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
+ call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
+ buftubebot=bordtubebot+tubebufthick
+ buftubetop=bordtubetop-tubebufthick
+ endif
+c if (shield_mode.gt.0) then
+c pi=3.141592d0
+C VSolvSphere the volume of solving sphere
+C print *,pi,"pi"
+C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
+C there will be no distinction between proline peptide group and normal peptide
+C group in case of shielding parameters
+c write (iout,*) "rpp(1,1)",rpp(1,1)," pi",pi
+c VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+c VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+c write (iout,*) "VSolvSphere",VSolvSphere,"VSolvSphere_div",
+c & VSolvSphere_div
+C long axis of side chain
+c do i=1,ntyp
+c long_r_sidechain(i)=vbldsc0(1,i)
+c short_r_sidechain(i)=sigma0(i)
+c enddo
+c buff_shield=1.0d0
+c endif
if (me.eq.king .or. .not.out1file )
& write (iout,*) "DISTCHAINMAX",distchainmax
do i=1,nrep
iremd_m_total=iremd_m_total+remd_m(i)
enddo
- write (iout,*) 'Total number of replicas ',iremd_m_total
+ write (iout,*) 'Total number of replicas ',iremd_m_total
+ endif
endif
- endif
+ if (adaptive) then
+ write (iout,*)
+ & "Adaptive (PMF-biased) umbrella sampling will be run"
+ call PMFread
+ endif
if(me.eq.king.or..not.out1file)
- & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
+ & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
return
end
c--------------------------------------------------------------------------
include 'COMMON.SETUP'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
+ include 'COMMON.FFIELD'
character*80 ucase
character*320 controlcard
ntime_split0=ntime_split
call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
ntime_split0=ntime_split
- call reada(controlcard,"R_CUT",r_cut,2.0d0)
- call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+c call reada(controlcard,"R_CUT",r_cut,2.0d0)
+c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
rest = index(controlcard,"REST").gt.0
tbf = index(controlcard,"TBF").gt.0
usampl = index(controlcard,"USAMPL").gt.0
+ scale_umb = index(controlcard,"SCALE_UMB").gt.0
+ adaptive = index(controlcard,"ADAPTIVE").gt.0
mdpdb = index(controlcard,"MDPDB").gt.0
call reada(controlcard,"T_BATH",t_bath,300.0d0)
call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
c if performing umbrella sampling, fragments constrained are read from the fragment file
nset=0
if(usampl) then
+ write (iout,*) "Umbrella sampling will be run"
+ if (scale_umb.and.adaptive) then
+ write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
+ write (iout,*) "Select one of those and re-run the job."
+ stop
+ endif
+ if (scale_umb) write (iout,*)
+ &"Umbrella-restraint force constants will be scaled by temperature"
call read_fragments
endif
enddo
write(iout,*) "angle constraints within ", nfrag_back,
& "backbone fragments."
+ if (loc_qlike) then
+ write(iout,*) "fragment, weights, q0:"
+ do i=1,nfrag_back
+ write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
+ & ifrag_back(2,i,iset),
+ & wfrag_back(1,i,iset),qin_back(1,i,iset),
+ & wfrag_back(2,i,iset),qin_back(2,i,iset),
+ & wfrag_back(3,i,iset),qin_back(3,i,iset)
+ enddo
+ else
write(iout,*) "fragment, weights:"
do i=1,nfrag_back
write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
& ifrag_back(2,i,iset),wfrag_back(1,i,iset),
& wfrag_back(2,i,iset),wfrag_back(3,i,iset)
enddo
+ endif
enddo
iset=mod(kolor,nset)+1
endif
include 'COMMON.BOUNDS'
include 'COMMON.MD'
include 'COMMON.SETUP'
+ include 'COMMON.SHIELD'
character*4 sequence(maxres)
integer rescode
double precision x(maxvar)
character*256 pdbfile
- character*320 weightcard
+ character*400 weightcard
character*80 weightcard_t,ucase
dimension itype_pdb(maxres)
common /pizda/ itype_pdb
logical seq_comp,fail
double precision energia(0:n_ene)
+ double precision secprob(3,maxdih_constr)
integer ilen
external ilen
C
C Body
C
C Read weights of the subsequent energy terms.
- 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
+ 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)
+ call reada(weightcard,'WSHIELD',wshield,1.0d0)
+ call reada(weightcard,'WLT',wliptran,0.0D0)
+ call reada(weightcard,'WTUBE',wtube,1.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
+ 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
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,
& 'General scaling factor of SC-p interactions:',scalscp
endif
r0_corr=cutoff_corr-delt_corr
- do i=1,20
+ do i=1,ntyp
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
+ wumb=1.0d0
call rescale_weights(t_bath)
if(me.eq.king.or..not.out1file)
& write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
call reada(weightcard,"V2SS",v2ss,7.61d0)
call reada(weightcard,"V3SS",v3ss,13.7d0)
call reada(weightcard,"EBR",ebr,-5.50D0)
+ call reada(weightcard,"ATRISS",atriss,0.301D0)
+ call reada(weightcard,"BTRISS",btriss,0.021D0)
+ call reada(weightcard,"CTRISS",ctriss,1.001D0)
+ call reada(weightcard,"DTRISS",dtriss,1.001D0)
+ write (iout,*) "ATRISS=", atriss
+ write (iout,*) "BTRISS=", btriss
+ write (iout,*) "CTRISS=", ctriss
+ write (iout,*) "DTRISS=", dtriss
+ dyn_ss=(index(weightcard,'DYN_SS').gt.0)
+ do i=1,maxres
+ dyn_ss_mask(i)=.false.
+ enddo
+ do i=1,maxres-1
+ do j=i+1,maxres
+ dyn_ssbond_ij(i,j)=1.0d300
+ enddo
+ enddo
+ call reada(weightcard,"HT",Ht,0.0D0)
+ if (dyn_ss) then
+ ss_depth=ebr/wsc-0.25*eps(1,1)
+ Ht=Ht/wsc-0.25*eps(1,1)
+ akcm=akcm*wstrain/wsc
+ akth=akth*wstrain/wsc
+ akct=akct*wstrain/wsc
+ v1ss=v1ss*wstrain/wsc
+ v2ss=v2ss*wstrain/wsc
+ v3ss=v3ss*wstrain/wsc
+ else
+ if (wstrain.ne.0.0) then
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ else
+ ss_depth=0.0
+ endif
+ endif
+ write (iout,*) "wshield,", wshield
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
-c print *,'indpdb=',indpdb,' pdbref=',pdbref
+ write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
+ write (iout,*)" HT",Ht
+ print *,'indpdb=',indpdb,' pdbref=',pdbref
endif
if (indpdb.gt.0 .or. pdbref) then
read(inp,'(a)') pdbfile
33 write (iout,'(a)') 'Error opening PDB file.'
stop
34 continue
-c print *,'Begin reading pdb data'
+c write (iout,*) 'Begin reading pdb data'
+c call flush(iout)
call readpdb
-c print *,'Finished reading pdb data'
+ nres_pdb=nres
+c write (iout,*) 'Finished reading pdb data'
+c call flush(iout)
if(me.eq.king.or..not.out1file)
& write (iout,'(a,i3,a,i3)')'nsup=',nsup,
& ' nstart_sup=',nstart_sup
maxsi=1000
do i=2,nres-1
iti=itype(i)
- if (iti.ne.10 .and. itype(i).ne.21) then
+ if (iti.ne.10 .and. itype(i).ne.ntyp1) then
nsi=0
fail=.true.
do while (fail.and.nsi.le.maxsi)
enddo
endif
endif
- if (indpdb.eq.0) then
+ if (indpdb.lt.2) then
C Read sequence if not taken from the pdb file.
read (inp,*) nres
c print *,'nres=',nres
vbld_inv(i)=vblinv
enddo
do i=2,nres-1
- vbld(i+nres)=dsc(itype(i))
- vbld_inv(i+nres)=dsc_inv(itype(i))
+ vbld(i+nres)=dsc(iabs(itype(i)))
+ vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
c write (iout,*) "i",i," itype",itype(i),
c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
enddo
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
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
#else
- if (itype(i).eq.21) then
+ if (itype(i).eq.ntyp1) then
#endif
itel(i)=0
#ifdef PROCOR
- else if (itype(i+1).ne.20) then
+ else if (iabs(itype(i+1)).ne.20) then
#else
- else if (itype(i).ne.20) then
+ else if (iabs(itype(i)).ne.20) then
#endif
itel(i)=1
else
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 (itype(1).eq.ntyp1) nnt=2
+ if (itype(nres).eq.ntyp1) nct=nct-1
+ if (indpdb.eq.1 .or. pdbref) then
if(me.eq.king.or..not.out1file)
& write (iout,'(a,i3)') 'nsup=',nsup
nstart_seq=nnt
if(me.eq.king.or..not.out1file)
& write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
& ' nstart_seq=',nstart_seq
+c AL 7/13/18 Compute quantities to locate chain
+ do i=nstart_sup,nstart_sup+nsup-2
+ call refsys(i,i+1,i+2,prod(1,1,i+2),prod(1,2,i+2),
+ & prod(1,3,i+2),fail)
+ if (fail) then
+ write (iout,*)
+ & "Failed to calculate reference system, residue",i
+ stop
+ endif
+ enddo
+ 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
+ write (iout,*) "ndih_constr",ndih_constr
+ if (ndih_constr.gt.0) then
+ raw_psipred=.false.
+C read (inp,*) ftors
+ read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
+ & i=1,ndih_constr)
+#ifdef MPI
+ if(me.eq.king.or..not.out1file)then
+#endif
+ write (iout,*)
+ & 'There are',ndih_constr,' restraints on gamma angles.'
+ do i=1,ndih_constr
+ write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
+ & ftors(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+ do i=1,ndih_constr
+ phi0(i)=deg2rad*phi0(i)
+ drange(i)=deg2rad*drange(i)
+ enddo
+C if(me.eq.king.or..not.out1file)
+C & 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
+#ifdef MPI
+ if (me.eq.king .or. .not.out1file) then
+#endif
+ write (iout,'(a)') 'Boundaries in gamma 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
+ else if (ndih_constr.lt.0) then
+ raw_psipred=.true.
+ call card_concat(weightcard)
+ call reada(weightcard,"PHIHEL",phihel,50.0D0)
+ call reada(weightcard,"PHIBET",phibet,180.0D0)
+ call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
+ call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
+ call reada(weightcard,"WDIHC",wdihc,0.591D0)
+ write (iout,*) "Weight of dihedral angle restraints",wdihc
+ read(inp,'(9x,3f7.3)')
+ & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
+ write (iout,*) "The secprob array"
+ do i=nnt,nct
+ write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
+ enddo
+ ndih_constr=0
+ do i=nnt+3,nct
+ if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
+ & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
+ ndih_constr=ndih_constr+1
+ idih_constr(ndih_constr)=i
+ sumv=0.0d0
+ do j=1,3
+ vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
+ sumv=sumv+vpsipred(j,ndih_constr)
+ enddo
+ do j=1,3
+ vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
+ enddo
+ phibound(1,ndih_constr)=phihel*deg2rad
+ phibound(2,ndih_constr)=phibet*deg2rad
+ sdihed(1,ndih_constr)=sigmahel*deg2rad
+ sdihed(2,ndih_constr)=sigmabet*deg2rad
+ endif
+ enddo
+#ifdef MPI
+ if(me.eq.king.or..not.out1file)then
+#endif
+ write (iout,*)
+ & 'There are',ndih_constr,
+ & ' bimodal restraints on gamma angles.'
+ do i=1,ndih_constr
+ write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
+ & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
+ & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
+ & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
+ & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
+ & (vpsipred(j,i),j=1,3)
+ enddo
+#ifdef MPI
+ endif
+#endif
+c Raw psipred input
endif
+C first setting the theta boundaries to 0 to pi
+C this mean that there is no energy penalty for any angle occuring this can be applied
+C for generate random conformation but is not implemented in this way
+C do i=1,nres
+C thetabound(1,i)=0
+C thetabound(2,i)=pi
+C enddo
+C begin reading theta constrains this is quartic constrains allowing to
+C have smooth second derivative
+ if (with_theta_constr) then
+C with_theta_constr is keyword allowing for occurance of theta constrains
+ read (inp,*) ntheta_constr
+C ntheta_constr is the number of theta constrains
+ if (ntheta_constr.gt.0) then
+C read (inp,*) ftors
+ read (inp,*) (itheta_constr(i),theta_constr0(i),
+ & theta_drange(i),for_thet_constr(i),
+ & i=1,ntheta_constr)
+C the above code reads from 1 to ntheta_constr
+C itheta_constr(i) residue i for which is theta_constr
+C theta_constr0 the global minimum value
+C theta_drange is range for which there is no energy penalty
+C for_thet_constr is the force constant for quartic energy penalty
+C E=k*x**4
+ if(me.eq.king.or..not.out1file)then
+ write (iout,*)
+ & 'There are',ntheta_constr,' constraints on phi angles.'
+ do i=1,ntheta_constr
+ write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
+ & theta_drange(i),
+ & for_thet_constr(i)
+ enddo
+ endif
+ do i=1,ntheta_constr
+ theta_constr0(i)=deg2rad*theta_constr0(i)
+ theta_drange(i)=deg2rad*theta_drange(i)
+ enddo
+C if(me.eq.king.or..not.out1file)
+C & write (iout,*) 'FTORS',ftors
+C do i=1,ntheta_constr
+C ii = itheta_constr(i)
+C thetabound(1,ii) = phi0(i)-drange(i)
+C thetabound(2,ii) = phi0(i)+drange(i)
+C enddo
+ endif ! ntheta_constr.gt.0
+ endif! with_theta_constr
+C
+C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
+C write (iout,*) "with_dihed_constr ",with_dihed_constr
c--- Zscore rms -------
if (nz_start.eq.0) nz_start=nnt
if (nz_end.eq.0 .and. nsup.gt.0) then
call MPI_Finalize(MPI_COMM_WORLD,IERROR)
stop 'Error reading reference structure'
#endif
- 39 call chainbuild
+ 39 call chainbuild_extconf
call setup_var
czscore call geom_to_var(nvar,coord_exp_zs(1,1))
nstart_sup=nnt
enddo
call contact(.true.,ncont_ref,icont_ref,co)
endif
-c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+ endif
+c print *, "A TU"
+ write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
call flush(iout)
if (constr_dist.gt.0) call read_dist_constr
write (iout,*) "After read_dist_constr nhpb",nhpb
+ if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
call hpb_partition
if(me.eq.king.or..not.out1file)
& write (iout,*) 'Contact order:',co
& restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
enddo
endif
- endif
- if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
+C endif
+ if (indpdb.lt.2 .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 (iranconf.eq.0 .and. indpdb.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
& ((c(l,k),l=1,3),k=1,nres),
& ((c(l,k+nres),l=1,3),k=nnt,nct)
write (iout,*) "Exit READ_CART"
- write (iout,'(8f10.5)')
- & ((c(l,k),l=1,3),k=1,nres),
- & ((c(l,k+nres),l=1,3),k=nnt,nct)
- call int_from_cart1(.true.)
- write (iout,*) "Finish INT_TO_CART"
+c write (iout,'(8f10.5)')
+c & ((c(l,k),l=1,3),k=1,nres),
+c & ((c(l,k+nres),l=1,3),k=nnt,nct)
+ call cartprint
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)
+c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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)
+c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ else
+ do j=1,3
+ dc(j,i+nres)=0.0d0
+c dc_norm(j,i+nres)=0.0d0
enddo
endif
enddo
- return
+ call int_from_cart1(.true.)
+ write (iout,*) "Finish INT_TO_CART"
+c write (iout,*) "DC"
+c do i=1,nres
+c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+c & (dc(j,i+nres),j=1,3)
+c enddo
+c call cartprint
+c call setup_var
+c return
else
call read_angles(inp,*36)
+ call chainbuild_extconf
endif
goto 37
36 write (iout,'(a)') 'Error reading angle file.'
enddo
do i=2,nres-1
omeg(i)=-120d0*deg2rad
+ if (itype(i).le.0) omeg(i)=-omeg(i)
enddo
+ call chainbuild_extconf
else
- if(me.eq.king.or..not.out1file)
- & write (iout,'(a)') 'Random-generated initial geometry.'
-
+ if (indpdb.eq.0) then
+ if (me.eq.king.or..not.out1file)
+ & write (iout,'(a)') 'Random-generated initial geometry.'
+ else
+ if (me.eq.king.or..not.out1file)write (iout,'(a)')
+ &'The geometry of missing residues in PDB file randomly generated.'
+ endif
#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
+ if (indpdb.eq.0) then
+ itmp=1
+ nrestemp=nres
+ call gen_rand_conf(itmp,nrestemp,*30)
+ goto 40
+ else
+ c = 0.0d0
+ do i=nstart_seq,nstart_seq+nsup
+ do j=1,3
+ c(j,i)=cref(j,i-nstart_seq+nnt,1)
+ c(j,i+nres)=cref(j,i-nstart_seq+nnt+nres_pdb,1)
+ enddo
+ enddo
+ write (iout,*)
+ & "Coordinates copied from the reference structure"
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
+ & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+ & (c(j,nres+ires),j=1,3)
+ enddo
+ if (nstart_seq.gt.nnt) then
+ nstartmp=nstart_seq-1
+ itmp=1
+ write (iout,*) "Building a random chain from residue",
+ & nstartmp," down to residue",itmp
+ call gen_rand_conf(nstartmp,itmp,*30)
+ goto 40
+ endif
+ itmp=nstart_seq+nsup+1
+ if (itmp.lt.nres) then
+ nrestemp=nres
+ write (iout,*)"Building a random chain from residue",
+ & itmp," up to residue",nrestemp
+ call gen_rand_conf(itmp,nrestemp,*30)
+ goto 40
+ endif
+ endif
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
40 continue
endif
#else
- 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 (*,*) 'Failed to generate random conformation',
- & ', itrial=',itrial
- enddo
- write (iout,'(a,i3,a)') 'Processor:',me,
- & ' error in generating random conformation.'
- write (*,'(a,i3,a)') 'Processor:',me,
+ write (*,'(a)')
& ' error in generating random conformation.'
stop
40 continue
#endif
+ write (iout,*) "Coordinates after building missing part"
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
+ & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+ & (c(j,nres+ires),j=1,3)
+ enddo
+ call int_from_cart1(.true.)
endif
elseif (modecalc.eq.4) then
read (inp,'(a)') intinname
write (iout,'(/a,i3,a)')
& 'The chain contains',ns,' disulfide-bridging cysteines.'
write (iout,'(20i4)') (iss(i),i=1,ns)
+ if (dyn_ss) then
+ write(iout,*)"Running with dynamic disulfide-bond formation"
+ else
write (iout,'(/a/)') 'Pre-formed links are:'
do i=1,nss
i1=ihpb(i)-nres
i2=jhpb(i)-nres
it1=itype(i1)
it2=itype(i2)
- if (me.eq.king.or..not.out1file)
- & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
+ write (iout,'(2a,i3,3a,i3,a,3f10.3)')
& restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
& ebr,forcon(i)
enddo
write (iout,'(a)')
+ endif
+ endif
+ if (ns.gt.0.and.dyn_ss) then
+ do i=nss+1,nhpb
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+ enddo
endif
+c --- test
+c call cartprint
+c write (iout,*) "DC"
+c do i=1,nres
+c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+c & (dc(j,i+nres),j=1,3)
+c enddo
if (i2ndstr.gt.0) call secstrp2dihc
c call geom_to_var(nvar,x)
c call etotal(energia(0))
& write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
& 'Processor',myrank,': end reading molecular data.'
#endif
+c print *,"A TU?"
return
end
c--------------------------------------------------------------------------
include 'COMMON.SETUP'
C Read bridging residues.
read (inp,*) ns,(iss(i),i=1,ns)
-c print *,'ns=',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 (iabs(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,
+ & 'Do you REALLY think that the residue ',
+ & restyp(itype(iss(i))),i,
& ' can form a disulfide bridge?!!!'
write (*,'(2a,i3,a)')
- & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
+ & 'Do you REALLY think that the residue ',
+ & restyp(itype(iss(i))),i,
& ' can form a disulfide bridge?!!!'
#ifdef MPI
call MPI_Finalize(MPI_COMM_WORLD,ierror)
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(fg_rank.eq.0)
+ & 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
enddo
enddo
do i=nnt,nct
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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)
nphi=nres-3
nvar=ntheta+nphi
nside=0
+ write (iout,*) "SETUP_VAR ialph"
do i=2,nres-1
- if (itype(i).ne.10 .and. itype(i).ne.21) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
nside=nside+1
ialph(i,1)=nvar+nside
ialph(nside,2)=i
else
nvar=nvar+2*nside
endif
-cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
+ write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
return
end
c----------------------------------------------------------------------------
character*3 out1file_text,ucase
character*3 ll
external ucase
+ tmpdir=""
+ out1file_text="YES"
c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
call getenv_loc("PREFIX",prefix)
pref_orig = prefix
open (ielep,file=elename,status='old',readonly,shared)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly,shared)
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',readonly,shared)
#elif (defined CRAY) || (defined AIX)
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
& action='read')
c print *,"Processor",myrank," opened file IBOND"
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old',action='read')
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old',action='read')
+#endif
c print *,"Processor",myrank," opened file ITHEP"
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old',action='read')
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',action='read')
+#endif
c print *,"Processor",myrank," opened file IROTAM"
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old',action='read')
c print *,"Processor",myrank," opened file IELEP"
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',action='read')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
c print *,"Processor",myrank," opened file ISIDEP"
c print *,"Processor",myrank," opened parameter files"
#elif (defined G77)
open (ibond,file=bondname,status='old')
call getenv_loc('THETPAR',thetname)
open (ithep,file=thetname,status='old')
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old')
+#endif
call getenv_loc('ROTPAR',rotname)
open (irotam,file=rotname,status='old')
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old')
+#endif
call getenv_loc('TORPAR',torname)
open (itorp,file=torname,status='old')
call getenv_loc('TORDPAR',tordname)
open (ielep,file=elename,status='old')
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old')
#else
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
& readonly)
open (itordp,file=tordname,status='old',readonly)
call getenv_loc('SCCORPAR',sccorname)
open (isccor,file=sccorname,status='old',readonly)
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ open (ithep_pdb,file=thetname_pdb,status='old',action='read')
+#endif
call getenv_loc('FOURIER',fouriername)
open (ifourier,file=fouriername,status='old',readonly)
call getenv_loc('ELEPAR',elename)
open (ielep,file=elename,status='old',readonly)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly)
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',action='read')
+#endif
+#endif
+#ifdef TUBE
+ call getenv_loc('TUBEPAR',tubename)
+#if defined(WINIFL) || defined(WINPGI)
+ open (itube,file=tubename,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open (itube,file=tubename,status='old',action='read')
+#elif (defined G77)
+ open (itube,file=tubename,status='old')
+#else
+ open (itube,file=tubename,status='old',readonly)
+#endif
#endif
#ifndef OLDSCP
C
mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
if (lentmp.gt.0)
- & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
& //'.stat')
rest2name=prefix(:ilen(prefix))//'.rst'
if(usampl) then
qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
endif
#endif
-#if defined(AIX) || defined(PGI)
- if (me.eq.king .or. .not. out1file)
- & open(iout,file=outname,status='unknown')
+#if defined(AIX) || defined(PGI) || defined(CRAY)
+ if (me.eq.king .or. .not. out1file) then
+ open(iout,file=outname,status='unknown')
+ else
+ open(iout,file="/dev/null",status="unknown")
+ endif
#ifdef DEBUG
if (fg_rank.gt.0) then
write (liczba,'(i3.3)') myrank/nfgtasks
& sccorname(:ilen(sccorname))
write (iout,*) "Bond & inertia constant file : ",
& bondname(:ilen(bondname))
- write (iout,*) "Bending parameter file : ",
+ write (iout,*) "Bending-torsion parameter file : ",
& thetname(:ilen(thetname))
write (iout,*) "Rotamer parameter file : ",
& rotname(:ilen(rotname))
+ write (iout,*) "Thetpdb parameter file : ",
+ & thetname_pdb(:ilen(thetname_pdb))
write (iout,*) "Threading database : ",
& patname(:ilen(patname))
if (lentmp.ne.0)
include 'COMMON.MD'
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
+ totTafm=totT
do i=1,2*nres
read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
enddo
include 'COMMON.MD'
include 'COMMON.CONTROL'
read(inp,*) nset,nfrag,npair,nfrag_back
+ loc_qlike=(nfrag_back.lt.0)
+ nfrag_back=iabs(nfrag_back)
if(me.eq.king.or..not.out1file)
& write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
- & " nfrag_back",nfrag_back
+ & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
do iset=1,nset
read(inp,*) mset(iset)
do i=1,nfrag
& write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
& ipair(2,i,iset), qinpair(i,iset)
enddo
+ if (loc_qlike) then
+ do i=1,nfrag_back
+ read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
+ & wfrag_back(2,i,iset),qin_back(2,i,iset),
+ & wfrag_back(3,i,iset),qin_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),qin_back(2,i,iset),
+ & wfrag_back(2,i,iset),qin_back(3,i,iset),
+ & wfrag_back(3,i,iset),qin_back(3,i,iset),
+ & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ enddo
+ else
do i=1,nfrag_back
read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
& wfrag_back(3,i,iset),
& 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
+ endif
enddo
return
end
+C---------------------------------------------------------------------------
+ subroutine read_afminp
+ 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'
+ character*320 afmcard
+ print *, "wchodze"
+ call card_concat(afmcard)
+ call readi(afmcard,"BEG",afmbeg,0)
+ call readi(afmcard,"END",afmend,0)
+ call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
+ call reada(afmcard,"VEL",velAFMconst,0.0d0)
+ print *,'FORCE=' ,forceAFMconst
+CCCC NOW PROPERTIES FOR AFM
+ distafminit=0.0d0
+ do i=1,3
+ distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
+ enddo
+ distafminit=dsqrt(distafminit)
+ print *,'initdist',distafminit
+ return
+ end
+
c-------------------------------------------------------------------------------
subroutine read_dist_constr
implicit real*8 (a-h,o-z)
integer ifrag_(2,100),ipair_(2,100)
double precision wfrag_(100),wpair_(100)
character*500 controlcard
-c write (iout,*) "Calling read_dist_constr"
+ logical normalize
+c print *, "WCHODZE"
+ 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 multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+ normalize = index(controlcard,"NORMALIZE").gt.0
c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
c write (iout,*) "IFRAG"
c do i=1,nfrag_
c 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
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
dhpb(nhpb)=ddjk
- forcon(nhpb)=wfrag_(i)
+ forcon(nhpb)=wfrag_(i)
else if (constr_dist.eq.2) then
if (ddjk.le.dist_cut) then
nhpb=nhpb+1
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#else
- write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
enddo
dhpb(nhpb)=dist(j,k)
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#else
- write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
enddo
enddo
endif
enddo
+c print *,ndist_
+ write (iout,*) "Distance restraints as read from input"
do i=1,ndist_
- read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
- if (forcon(nhpb+1).gt.0.0d0) then
+ if (constr_dist.eq.11) then
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+ & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+ if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
+ nhpb=nhpb+1
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king)
+ & write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb)
+#else
+ write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb)
+#endif
+ if (ibecarb(i).gt.0) then
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ endif
+ else
+C print *,"in else"
+ 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
- dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ 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
#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)
+ & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
#else
- write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
- & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
#endif
endif
+C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+C if (forcon(nhpb+1).gt.0.0d0) then
+C nhpb=nhpb+1
+C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
enddo
+ if (constr_dist.eq.11 .and. normalize) then
+ fordepthmax=fordepth(1)
+ do i=2,nhpb
+ if (fordepth(i).gt.fordepthmax) fordepthmax=fordepth(i)
+ enddo
+ do i=1,nhpb
+ fordepth(i)=fordepth(i)/fordepthmax
+ enddo
+ write (iout,'(/a/4a5,2a8,2a10)')
+ & "Normalized Lorenzian-like distance restraints",
+ & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V"
+ do i=1,nhpb
+ write (iout,'(4i5,2f8.2,2f10.5)')i,ihpb(i),jhpb(i),ibecarb(i),
+ & dhpb(i),dhpb1(i),forcon(i),fordepth(i)
+ enddo
+ endif
+ write (iout,*)
call flush(iout)
return
end
#endif
if (OKRandom) then
c r1 = prng_next(me)
- r1=ran_number(0.0D0,1.0D0)
+ r1=ran_number(0.0D0,1.0D0)
if(me.eq.king)
& write (iout,*) 'ran_num',r1
if (r1.lt.0.0d0) OKRandom=.false.