#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,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
- enddo
+ if (nhpb.gt.nss) then
+ write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
+ & "The following",nhpb-nss,
+ & " distance restraints have been imposed:",
+ & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
+ & " score"," type"
+ do i=nss+1,nhpb
+ write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
+ & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
+ & irestr_type(i)
+ enddo
+ endif
+ if (npeak.gt.0) then
+ write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
+ & "The following",npeak,
+ & " NMR peak restraints have been imposed:",
+ & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
+ & " score"," type"," ipeak"
+ do i=1,npeak
+ do j=ipeak(1,i),ipeak(2,i)
+ write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
+ & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
+ & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
+ enddo
+ enddo
+ write (iout,*) "The ipeak array"
+ do i=1,npeak
+ write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
+ enddo
+ endif
#ifdef MPI
endif
#endif
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)
+ write (iout,*) "constr_dist",constr_dist
+ call readi(controlcard,'NSAXS',nsaxs,0)
+ call readi(controlcard,'SAXS_MODE',saxs_mode,0)
+ call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
+ call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
+ write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+ & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
+ read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
call readi(controlcard,'SYM',symetr,1)
call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
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
+c write(iout,*) "bordliptop=",bordliptop
+c write(iout,*) "bordlipbot=",bordlipbot
+c write(iout,*) "bufliptop=",bufliptop
+c write(iout,*) "buflipbot=",buflipbot
if (me.eq.king .or. .not.out1file )
call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
call reada(weightcard,'TEMP0',temp0,300.0d0)
call reada(weightcard,'WLT',wliptran,0.0D0)
+ call reada(weightcard,'WSAXS',wsaxs,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)
weights(17)=wbond
weights(18)=scal14
weights(21)=wsccor
+ weights(25)=wsaxs
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,
call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
nsi=nsi+1
enddo
+c AL 7/10/16
+c Calculalte only the coordinates of the current sidechain; no need to rebuild
+c whole chain
+ call locate_side_chain(i)
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 call chainbuild
endif
endif
if (indpdb.eq.0) then
phibound(2,i)=pi
enddo
read (inp,*) ndih_constr
+ write (iout,*) "ndish_constr",ndih_constr
if (ndih_constr.gt.0) then
read (inp,*) ftors
read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
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
& restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
enddo
endif
- endif
-
-
+ write (iout,*) "calling read_saxs_consrtr",nsaxs
+ if (nsaxs.gt.0) call read_saxs_constr
+
if (constr_homology.gt.0) then
call read_constr_homology
if (indpdb.gt.0 .or. pdbref) then
if (nhpb.gt.0) call hpb_partition
+ if (peak.gt.0) call NMRpeak_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
return
else
+ write (iout,*) "Calling read_ang"
call read_angles(inp,*36)
+ write (iout,*) "Calling chainbuild"
+ call chainbuild_extconf
endif
goto 37
36 write (iout,'(a)') 'Error reading angle file.'
omeg(i)=-120d0*deg2rad
if (itype(i).le.0) omeg(i)=-omeg(i)
enddo
+c from old chainbuild
+C
+C Define the origin and orientation of the coordinate system and locate the
+C first three CA's and SC(2).
+C
+ call orig_frame
+*
+* Build the alpha-carbon chain.
+*
+ do i=4,nres
+ call locate_next_res(i)
+ enddo
+C
+C First and last SC must coincide with the corresponding CA.
+C
+ do j=1,3
+ dc(j,nres+1)=0.0D0
+ dc_norm(j,nres+1)=0.0D0
+ dc(j,nres+nres)=0.0D0
+ dc_norm(j,nres+nres)=0.0D0
+ c(j,nres+1)=c(j,1)
+ c(j,nres+nres)=c(j,nres)
+ enddo
+C
+C Define the origin and orientation of the coordinate system and locate the
+C first three CA's and SC(2).
+C
+ call orig_frame
+*
+* Build the alpha-carbon chain.
+*
+ do i=4,nres
+ call locate_next_res(i)
+ enddo
+C
+C First and last SC must coincide with the corresponding CA.
+C
+ do j=1,3
+ dc(j,nres+1)=0.0D0
+ dc_norm(j,nres+1)=0.0D0
+ dc(j,nres+nres)=0.0D0
+ dc_norm(j,nres+nres)=0.0D0
+ c(j,nres+1)=c(j,1)
+ c(j,nres+nres)=c(j,nres)
+ enddo
+
+c
else
if(me.eq.king.or..not.out1file)
& write (iout,'(a)') 'Random-generated initial geometry.'
include 'DIMENSIONS'
include 'COMMON.MINIM'
include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SETUP'
character*80 ucase
character*320 minimcard
call card_concat(minimcard)
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)
+#ifdef MPI
+ if (.not. out1file .or. me.eq.king) then
+#endif
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
+#ifdef MPI
+ endif
+#endif
return
end
c----------------------------------------------------------------------------
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)
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')
+ open (iliptranpar,file=liptranname,status='old',readonly)
#ifndef CRYST_SC
call getenv_loc('ROTPARPDB',rotname_pdb)
open (irotam_pdb,file=rotname_pdb,status='old',action='read')
end
c-------------------------------------------------------------------------------
+ subroutine read_saxs_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'
+ double precision cm(3)
+c read(inp,*) nsaxs
+ write (iout,*) "Calling read_saxs nsaxs",nsaxs
+ call flush(iout)
+ if (saxs_mode.eq.0) then
+c SAXS distance distribution
+ do i=1,nsaxs
+ read(inp,*) distsaxs(i),Psaxs(i)
+ enddo
+ Cnorm = 0.0d0
+ do i=1,nsaxs
+ Cnorm = Cnorm + Psaxs(i)
+ enddo
+ write (iout,*) "Cnorm",Cnorm
+ do i=1,nsaxs
+ Psaxs(i)=Psaxs(i)/Cnorm
+ enddo
+ write (iout,*) "Normalized distance distribution from SAXS"
+ do i=1,nsaxs
+ write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
+ enddo
+ Wsaxs0=0.0d0
+ do i=1,nsaxs
+ Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
+ enddo
+ write (iout,*) "Wsaxs0",Wsaxs0
+ else
+c SAXS "spheres".
+ do i=1,nsaxs
+ read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
+ enddo
+ do j=1,3
+ cm(j)=0.0d0
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ cm(j)=cm(j)+Csaxs(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/nsaxs
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ Csaxs(j,i)=Csaxs(j,i)-cm(j)
+ enddo
+ enddo
+ write (iout,*) "SAXS sphere coordinates"
+ do i=1,nsaxs
+ write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
+ enddo
+ endif
+ return
+ end
+c-------------------------------------------------------------------------------
subroutine read_dist_constr
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
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"
+ integer ifrag_(2,100),ipair_(2,1000)
+ double precision wfrag_(100),wpair_(1000)
+ character*5000 controlcard
+ logical normalize,next
+ integer restr_type
+ double precision xlink(4,0:4) /
+c a b c sigma
+ & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
+ & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
+ & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
+ & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
+ & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
+c print *, "WCHODZE"
+ write (iout,*) "Calling read_dist_constr"
c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
c call flush(iout)
+ next=.true.
+
+ npeak=0
+ ipeak=0
+ nhpb_peak=0
+
+ DO WHILE (next)
+
call card_concat(controlcard)
+ next = index(controlcard,"NEXT").gt.0
+ call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
+ write (iout,*) "restr_type",restr_type
+ call readi(controlcard,"NFRAG",nfrag_,0)
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)
+ if (restr_type.eq.10)
+ & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
+ if (restr_type.eq.12)
+ & call reada(controlcard,'SCAL_PEAK',scal_peak,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)
+ normalize = index(controlcard,"NORMALIZE").gt.0
+ write (iout,*) "WBOLTZD",wboltzd
+ write (iout,*) "SCAL_PEAK",scal_peak
+ write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+ write (iout,*) "IFRAG"
+ do i=1,nfrag_
+ write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ enddo
+ write (iout,*) "IPAIR"
+ do i=1,npair_
+ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+ enddo
+ if (nfrag_.gt.0) write (iout,*)
+ & "Distance restraints as generated from reference structure"
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
+c call flush(iout)
+ if (wfrag_(i).eq.0.0d0) cycle
do j=ifrag_(1,i),ifrag_(2,i)-1
do k=j+1,ifrag_(2,i)
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
+ if (restr_type.eq.1) then
+ nhpb=nhpb+1
+ irestr_type(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
+ irestr_type(nhpb)=1
ihpb(nhpb)=j
jhpb(nhpb)=k
dhpb(nhpb)=ddjk
forcon(nhpb)=wfrag_(i)
endif
- else
+ else if (restr_type.eq.3) then
nhpb=nhpb+1
+ irestr_type(nhpb)=1
ihpb(nhpb)=j
jhpb(nhpb)=k
dhpb(nhpb)=ddjk
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
enddo
- endif
enddo
do i=1,npair_
- if (wpair_(i).gt.0.0d0) then
+ if (wpair_(i).eq.0.0d0) cycle
ii = ipair_(1,i)
jj = ipair_(2,i)
if (ii.gt.jj) then
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)
+ ddjk=dist(j,k)
+ if (restr_type.eq.1) then
+ nhpb=nhpb+1
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wpair_(i)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wpair_(i)
+ endif
+ else if (restr_type.eq.3) then
+ nhpb=nhpb+1
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
#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,1pe12.2)') "+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,1pe12.2)') "+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 (restr_type.eq.12) then
+ read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
+ & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
+ & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
+ & fordepth_peak(nhpb_peak+1),npeak
+c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
+c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
+c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
+c & fordepth_peak(nhpb_peak+1),npeak
+ if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
+ & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
+ nhpb_peak=nhpb_peak+1
+ irestr_type_peak(nhpb_peak)=12
+ if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
+ ipeak(2,npeak)=i
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king)
+ & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
+ & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
+ & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
+ & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
+#else
+ write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
+ & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
+ & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
+ & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
+#endif
+ if (ibecarb_peak(nhpb_peak).gt.0) then
+ ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
+ jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
+ endif
+ else if (restr_type.eq.11) then
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
+ & dhpb1(nhpb+1),ibecarb(nhpb+1),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
+ irestr_type(nhpb)=11
+#ifdef MPI
+ if (.not.out1file .or. me.eq.king)
+ & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
+#else
+ write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
+#endif
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
+ endif
+ else if (restr_type.eq.10) then
+c Cross-lonk Markov-like potential
+ call card_concat(controlcard)
+ call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
+ call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
+ ibecarb(nhpb+1)=0
+ if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
+ if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
+ if (index(controlcard,"ZL").gt.0) then
+ link_type=1
+ else if (index(controlcard,"ADH").gt.0) then
+ link_type=2
+ else if (index(controlcard,"PDH").gt.0) then
+ link_type=3
+ else if (index(controlcard,"DSS").gt.0) then
+ link_type=4
+ else
+ link_type=0
+ endif
+ call reada(controlcard,"AXLINK",dhpb(nhpb+1),
+ & xlink(1,link_type))
+ call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
+ & xlink(2,link_type))
+ call reada(controlcard,"CXLINK",fordepth(nhpb+1),
+ & xlink(3,link_type))
+ call reada(controlcard,"SIGMA",forcon(nhpb+1),
+ & xlink(4,link_type))
+ call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
+c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
+c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
+ if (forcon(nhpb+1).le.0.0d0 .or.
+ & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
nhpb=nhpb+1
- dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ irestr_type(nhpb)=10
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
+ 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,2f8.2,3f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
+ & irestr_type(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,2f8.2,3f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
+ & irestr_type(nhpb)
+#endif
+ else
+C print *,"in else"
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
+ & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ if (dhpb1(nhpb).eq.0.0d0) then
+ irestr_type(nhpb)=1
+ else
+ irestr_type(nhpb)=2
+ endif
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+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,4i5,f8.2,f10.1)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
+#else
+ 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
- call flush(iout)
+
+ ENDDO ! next
+
+ fordepthmax=0.0d0
+ if (normalize) then
+ do i=nss+1,nhpb
+ if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
+ & fordepthmax=fordepth(i)
+ enddo
+ do i=nss+1,nhpb
+ if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
+ enddo
+ endif
return
end
c-------------------------------------------------------------------------------
logical lprn /.true./
integer ilen
external ilen
+ logical liiflag
c
c FP - Nov. 2014 Temporary specifications for new vars
c
call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
+ dist1cut=(index(controlcard,'DIST1CUT').gt.0)
call readi(controlcard,"HOMOL_NSET",homol_nset,1)
read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
if(.not.read2sigma.and.start_from_model) then
- write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
+ & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
start_from_model=.false.
endif
- if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
+ if(start_from_model .and. (me.eq.king .or. .not. out1file))
+ & write(iout,*) 'START_FROM_MODELS is ON'
if(start_from_model .and. rest) then
- write(iout,*) 'START_FROM_MODELS is OFF'
- write(iout,*) 'remove restart keyword from input'
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) 'START_FROM_MODELS is OFF'
+ write(iout,*) 'remove restart keyword from input'
+ endif
endif
if (homol_nset.gt.1)then
call card_concat(controlcard)
lim_odl=0
lim_dih=0
c
- write(iout,*) 'nnt=',nnt,'nct=',nct
+c write(iout,*) 'nnt=',nnt,'nct=',nct
c
do i = nnt,nct
do k=1,constr_homology
enddo
enddo
+ if (read_homol_frag) then
+ call read_klapaucjusz
+ else
+
do k=1,constr_homology
read(inp,'(a)') pdbfile
-c Next stament causes error upon compilation (?)
-c if(me.eq.king.or. .not. out1file)
-c write (iout,'(2a)') 'PDB data will be read from file ',
-c & pdbfile(:ilen(pdbfile))
- write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
+ if(me.eq.king .or. .not. out1file)
+ & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
& pdbfile(:ilen(pdbfile))
open(ipdbin,file=pdbfile,status='old',err=33)
goto 34
rescore(k,i_tmp)=rescore_tmp
rescore2(k,i_tmp)=rescore2_tmp
rescore3(k,i_tmp)=rescore3_tmp
- write(iout,'(a7,i5,3f10.5,i5)') "rescore",
+ if (.not. out1file .or. me.eq.king)
+ & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
& i_tmp,rescore2_tmp,rescore_tmp,
& rescore3_tmp,idomain_tmp
else
distal=dsqrt(x12*x12+y12*y12+z12*z12)
c write (iout,*) k,i,j,distal,dist2_cut
+ if (dist1cut .and. k.gt.1) then
+ ii=ii+1
+ if (l_homo(1,ii)) then
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+ odl(k,ii)=distal
+ sigma_odl(k,ii)=sigma_odl(1,ii)
+ else
+ l_homo(k,ii)=.false.
+ endif
+ else
if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
& .and. distal.le.dist2_cut ) then
ii=ii+1
l_homo(k,ii)=.false.
endif
+ endif
enddo
enddo
lim_odl=ii
c
if (waga_dist.ne.0.0d0) then
ii=0
+ liiflag=.true.
do i=nnt,nct-2
do j=i+2,nct
ii=ii+1
- if (ii_in_use(ii).eq.0) then
- do ki=ii,lim_odl-1
- ires_homo(ki)=ires_homo(ki+1)
- jres_homo(ki)=jres_homo(ki+1)
- ii_in_use(ki)=ii_in_use(ki+1)
+ if (ii_in_use(ii).eq.0.and.liiflag) then
+ liiflag=.false.
+ iistart=ii
+ endif
+ if (ii_in_use(ii).ne.0.and..not.liiflag.or.
+ & .not.liiflag.and.ii.eq.lim_odl) then
+ if (ii.eq.lim_odl) then
+ iishift=ii-iistart+1
+ else
+ iishift=ii-iistart
+ endif
+ liiflag=.true.
+ do ki=iistart,lim_odl-iishift
+ ires_homo(ki)=ires_homo(ki+iishift)
+ jres_homo(ki)=jres_homo(ki+iishift)
+ ii_in_use(ki)=ii_in_use(ki+iishift)
do k=1,constr_homology
- odl(k,ki)=odl(k,ki+1)
- sigma_odl(k,ki)=sigma_odl(k,ki+1)
- l_homo(k,ki)=l_homo(k,ki+1)
+ odl(k,ki)=odl(k,ki+iishift)
+ sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
+ l_homo(k,ki)=l_homo(k,ki+iishift)
enddo
enddo
- ii=ii-1
- lim_odl=lim_odl-1
+ ii=ii-iishift
+ lim_odl=lim_odl-iishift
endif
enddo
enddo
endif
+
+ endif ! .not. klapaucjusz
+
if (constr_homology.gt.0) call homology_partition
if (constr_homology.gt.0) call init_int_table
c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
#endif
return
end
+
+c -----------------------------------------------------------------
+ subroutine read_klapaucjusz
+
+ 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'
+ include 'COMMON.NAMES'
+ character*256 fragfile
+ integer ninclust(maxclust),inclust(max_template,maxclust),
+ & nresclust(maxclust),iresclust(maxres,maxclust)
+
+ character*2 kic2
+ character*24 model_ki_dist, model_ki_angle
+ character*500 controlcard
+ integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
+ logical lprn /.true./
+ integer ilen
+ external ilen
+ logical liiflag
+c
+c
+ double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+ double precision, dimension (max_template,maxres) :: rescore
+ double precision, dimension (max_template,maxres) :: rescore2
+ character*24 pdbfile,tpl_k_rescore
+
+c
+c For new homol impl
+c
+ include 'COMMON.VAR'
+c
+ call getenv("FRAGFILE",fragfile)
+ open(ientin,file=fragfile,status="old",err=10)
+ read(ientin,*) constr_homology,nclust
+ l_homo = .false.
+ sigma_theta=0.0
+ sigma_d=0.0
+ sigma_dih=0.0
+c Read pdb files
+ do k=1,constr_homology
+ read(ientin,'(a)') pdbfile
+ write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ stop
+ 34 continue
+ unres_pdb=.false.
+ call readpdb_template(k)
+ do i=1,nres
+ rescore(k,i)=0.2d0
+ rescore2(k,i)=1.0d0
+ enddo
+ enddo
+c Read clusters
+ do i=1,nclust
+ read(ientin,*) ninclust(i),nresclust(i)
+ read(ientin,*) (inclust(k,i),k=1,ninclust(i))
+ read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
+ enddo
+c
+c Loop over clusters
+c
+ do l=1,nclust
+ do ll = 1,ninclust(l)
+
+ k = inclust(ll,l)
+ do i=1,nres
+ idomain(k,i)=0
+ enddo
+ do i=1,nresclust(l)
+ if (nnt.gt.1) then
+ idomain(k,iresclust(i,l)+1) = 1
+ else
+ idomain(k,iresclust(i,l)) = 1
+ endif
+ enddo
+c
+c Distance restraints
+c
+c ... --> odl(k,ii)
+C Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=chomo(j,i,k)
+c write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+
+ x12=c(1,i)-c(1,j)
+ y12=c(2,i)-c(2,j)
+ z12=c(3,i)-c(3,j)
+ distal=dsqrt(x12*x12+y12*y12+z12*z12)
+c write (iout,*) k,i,j,distal,dist2_cut
+
+ if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+ & .and. distal.le.dist2_cut ) then
+
+ ii=ii+1
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+
+c write (iout,*) "k",k
+c write (iout,*) "i",i," j",j," constr_homology",
+c & constr_homology
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+ odl(k,ii)=distal
+ if (read2sigma) then
+ sigma_odl(k,ii)=0
+ do ik=i,j
+ sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+ enddo
+ sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+ if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
+ & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+ else
+ if (odl(k,ii).le.dist_cut) then
+ sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
+ else
+#ifdef OLDSIGMA
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+#else
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+#endif
+ endif
+ endif
+ sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
+ else
+ ii=ii+1
+c l_homo(k,ii)=.false.
+ endif
+ enddo
+ enddo
+ lim_odl=ii
+ endif
+c
+c Theta, dihedral and SC retraints
+c
+ if (waga_angle.gt.0.0d0) then
+ do i = nnt+3,nct
+ if (idomain(k,i).eq.0) then
+c sigma_dih(k,i)=0.0
+ cycle
+ endif
+ dih(k,i)=phiref(i)
+ sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)+rescore(k,i-3))/4.0
+c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
+c & " sigma_dihed",sigma_dih(k,i)
+ if (sigma_dih(k,i).ne.0)
+ & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+ enddo
+ lim_dih=nct-nnt-2
+ endif
+
+ if (waga_theta.gt.0.0d0) then
+ do i = nnt+2,nct
+ if (idomain(k,i).eq.0) then
+c sigma_theta(k,i)=0.0
+ cycle
+ endif
+ thetatpl(k,i)=thetaref(i)
+ sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2))/3.0
+ if (sigma_theta(k,i).ne.0)
+ & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+ enddo
+ endif
+
+ if (waga_d.gt.0.0d0) then
+ do i = nnt,nct
+ if (itype(i).eq.10) cycle
+ if (idomain(k,i).eq.0 ) then
+c sigma_d(k,i)=0.0
+ cycle
+ endif
+ xxtpl(k,i)=xxref(i)
+ yytpl(k,i)=yyref(i)
+ zztpl(k,i)=zzref(i)
+ sigma_d(k,i)=rescore(k,i)
+ if (sigma_d(k,i).ne.0)
+ & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+ if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
+ enddo
+ endif
+ enddo ! l
+ enddo ! ll
+c
+c remove distance restraints not used in any model from the list
+c shift data in all arrays
+c
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ liiflag=.true.
+ do i=nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ if (ii_in_use(ii).eq.0.and.liiflag) then
+ liiflag=.false.
+ iistart=ii
+ endif
+ if (ii_in_use(ii).ne.0.and..not.liiflag.or.
+ & .not.liiflag.and.ii.eq.lim_odl) then
+ if (ii.eq.lim_odl) then
+ iishift=ii-iistart+1
+ else
+ iishift=ii-iistart
+ endif
+ liiflag=.true.
+ do ki=iistart,lim_odl-iishift
+ ires_homo(ki)=ires_homo(ki+iishift)
+ jres_homo(ki)=jres_homo(ki+iishift)
+ ii_in_use(ki)=ii_in_use(ki+iishift)
+ do k=1,constr_homology
+ odl(k,ki)=odl(k,ki+iishift)
+ sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
+ l_homo(k,ki)=l_homo(k,ki+iishift)
+ enddo
+ enddo
+ ii=ii-iishift
+ lim_odl=lim_odl-iishift
+ endif
+ enddo
+ enddo
+ endif
+
+ return
+ 10 stop "Error infragment file"
+ end
+c----------------------------------------------------------------------
projects / unres.git / blobdiff