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)
+ call readi(controlcard,'CONSTR_HOMOL',constr_homology,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(weightcard,'WTORD',wtor_d,1.0D0)
call reada(weightcard,'WANG',wang,1.0D0)
call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+ call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
+ call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
+ call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
+ call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
call reada(weightcard,'SCAL14',scal14,0.4D0)
call reada(weightcard,'SCALSCP',scalscp,1.0d0)
call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
weights(18)=scal14
weights(21)=wsccor
endif
+ weights(25)=wdfa_dist
+ weights(26)=wdfa_tor
+ weights(27)=wdfa_nei
+ weights(28)=wdfa_beta
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
+ & wturn4,wturn6,
+ & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
10 format (/'Energy-term weights (unscaled):'//
& 'WSCC= ',f10.6,' (SC-SC)'/
& 'WSCP= ',f10.6,' (SC-p)'/
& 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
& 'WTURN3= ',f10.6,' (turns, 3rd order)'/
& 'WTURN4= ',f10.6,' (turns, 4th order)'/
- & 'WTURN6= ',f10.6,' (turns, 6th order)')
+ & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+ & 'WDFA_D= ',f10.6,' (DFA, distance)' /
+ & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
+ & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
+ & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
if(me.eq.king.or..not.out1file)then
if (wcorr4.gt.0.0d0) then
write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
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
+ & wturn4,wturn6,
+ & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
22 format (/'Energy-term weights (scaled):'//
& 'WSCC= ',f10.6,' (SC-SC)'/
& 'WSCP= ',f10.6,' (SC-p)'/
& 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
& 'WTURN3= ',f10.6,' (turns, 3rd order)'/
& 'WTURN4= ',f10.6,' (turns, 4th order)'/
- & 'WTURN6= ',f10.6,' (turns, 6th order)')
+ & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+ & 'WDFA_D= ',f10.6,' (DFA, distance)' /
+ & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
+ & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
+ & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
if(me.eq.king.or..not.out1file)
& write (iout,*) "Reference temperature for weights calculation:",
& temp0
call reada(weightcard,"V2SS",v2ss,7.61d0)
call reada(weightcard,"V3SS",v3ss,13.7d0)
call reada(weightcard,"EBR",ebr,-5.50D0)
+ 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
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ endif
+
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
+ 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
call contact(.false.,ncont_ref,icont_ref,co)
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
& i,' after ',nsi,' trials'
endif
enddo
+C 10/03/12 Adam: Recalculate coordinates with new side chain positions
+ call chainbuild
endif
+C Following 2 lines for diagnostics; comment out if not needed
+c write (iout,*) "After sideadd"
+c call intout
endif
if (indpdb.eq.0) then
C Read sequence if not taken from the pdb file.
vbld_inv(i)=vblinv
enddo
do i=2,nres-1
- vbld(i+nres)=dsc(iabs(itype(i)))
- vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
+ 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
itel(i)=0
#ifdef PROCOR
- else if (iabs(itype(i+1)).ne.20) then
+ else if (itype(i+1).ne.20) then
#else
- else if (iabs(itype(i)).ne.20) then
+ else if (itype(i).ne.20) then
#endif
itel(i)=1
else
cd print *,'NNT=',NNT,' NCT=',NCT
if (itype(1).eq.21) nnt=2
if (itype(nres).eq.21) nct=nct-1
+
+C Bartek:READ init_vars
+C Initialize variables!
+C Juyong:READ read_info
+C READ fragment information!!
+C both routines should be in dfa.F file!!
+
+ if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
+ & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
+ call init_dfa_vars
+ print*, 'init_dfa_vars finished!'
+ call read_dfa_info
+ print*, 'read_dfa_info finished!'
+ endif
+C
if (pdbref) then
if(me.eq.king.or..not.out1file)
& write (iout,'(a,i3)') 'nsup=',nsup
c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
if (constr_dist.gt.0) then
call read_dist_constr
- call hpb_partition
endif
+
+
+ if (constr_homology.gt.0) then
+ call read_constr_homology
+ else
+ homol_nset=0
+ 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
enddo
do i=2,nres-1
omeg(i)=-120d0*deg2rad
- if (itype(i).le.0) omeg(i)=-omeg(i)
enddo
else
if(me.eq.king.or..not.out1file)
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
if (i2ndstr.gt.0) call secstrp2dihc
c call geom_to_var(nvar,x)
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,
+ & '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
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.MD'
+ include 'COMMON.CONTROL'
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
do i=1,2*nres
do i=1,2*nres
read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
read (irest2,*) iset
endif
close(irest2)
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
+c write (iout,*) "j",j," k",k
ddjk=dist(j,k)
if (constr_dist.eq.1) then
nhpb=nhpb+1
return
end
c-------------------------------------------------------------------------------
+
+ subroutine read_constr_homology
+
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MD'
+ include 'COMMON.GEO'
+ include 'COMMON.INTERACT'
+ double precision odl_temp,sigma_odl_temp
+ common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+ & sigma_odl_temp(maxres,maxres,max_template)
+ character*2 kic2
+ character*24 model_ki_dist, model_ki_angle
+ character*500 controlcard
+ character*3200 controlcard1
+ integer ki, i, j, k, l
+ logical lprn /.true./
+
+ call card_concat(controlcard)
+ call readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ if (homol_nset.gt.1)then
+ call card_concat(controlcard)
+ read(controlcard,*) (waga_dist(i),i=1,homol_nset)
+ call card_concat(controlcard)
+ read(controlcard,*) (waga_angle(i),i=1,homol_nset)
+ write(iout,*) "iset distance_weight angle_weight"
+ do i=1,homol_nset
+ write(iout,*) i,waga_dist(i),waga_angle(i)
+ enddo
+ else
+ iset=1
+ call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
+ call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
+ endif
+
+ write (iout,*) "nnt",nnt," nct",nct
+ call flush(iout)
+ lim_odl=0
+ lim_dih=0
+ do i=1,nres
+ do j=i+2,nres
+ do ki=1,constr_homology
+ sigma_odl_temp(i,j,ki)=0.0d0
+ odl_temp(i,j,ki)=0.0d0
+ enddo
+ enddo
+ enddo
+ do i=1,nres-3
+ do ki=1,constr_homology
+ dih(ki,i)=0.0d0
+ sigma_dih(ki,i)=0.0d0
+ enddo
+ enddo
+ do ki=1,constr_homology
+ write(kic2,'(i2)') ki
+ if (ki.le.9) kic2="0"//kic2(2:2)
+
+ model_ki_dist="model"//kic2//".dist"
+ model_ki_angle="model"//kic2//".angle"
+ open (ientin,file=model_ki_dist,status='old')
+ do irec=1,maxdim !petla do czytania wiezow na odleglosc
+ read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
+ & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
+ odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
+ sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
+ & sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
+ enddo
+ 1401 continue
+ close (ientin)
+ open (ientin,file=model_ki_angle,status='old')
+ do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
+ read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
+ & sigma_dih(ki,i+nnt-1)
+ if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
+ sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
+ enddo
+ 1402 continue
+ close (ientin)
+ enddo
+ ii=0
+ write (iout,*) "nnt",nnt," nct",nct
+ do i=nnt,nct-2
+ do j=i+2,nct
+ ki=1
+c write (iout,*) "i",i," j",j," constr_homology",constr_homology
+ do while (ki.le.constr_homology .and.
+ & sigma_odl_temp(i,j,ki).le.0.0d0)
+c write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
+ ki=ki+1
+ enddo
+c write (iout,*) "ki",ki
+ if (ki.gt.constr_homology) cycle
+ ii=ii+1
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+ do ki=1,constr_homology
+ odl(ki,ii)=odl_temp(i,j,ki)
+ sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
+ enddo
+ enddo
+ enddo
+ lim_odl=ii
+ if (constr_homology.gt.0) call homology_partition
+c Print restraints
+ if (.not.lprn) return
+ write (iout,*) "Distance restraints from templates"
+ do ii=1,lim_odl
+ write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
+ & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
+ enddo
+ write (iout,*) "Dihedral angle restraints from templates"
+ do i=nnt,lim_dih
+ write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+ & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
+ enddo
+c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
+c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
+
+
+ return
+ end
+c----------------------------------------------------------------------
+
#ifdef WINIFL
subroutine flush(iu)
return
return
end
#endif
+
c------------------------------------------------------------------------------
subroutine copy_to_tmp(source)
include "DIMENSIONS"
projects / unres.git / blobdiff