subroutine readrtns
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.SBRIDGE'
include 'COMMON.IOUNITS'
include 'COMMON.SPLITELE'
+ integer i,j
logical file_exist
C Read job setup parameters
call read_control
C
C Read contorl data
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MP
include 'mpif.h'
logical OKRandom, prng_restart
- real*8 r1
+ double precision r1
#endif
include 'COMMON.IOUNITS'
include 'COMMON.TIME1'
include 'COMMON.THREAD'
include 'COMMON.SBRIDGE'
include 'COMMON.CONTROL'
+ include 'COMMON.SAXS'
include 'COMMON.MCM'
include 'COMMON.MAP'
include 'COMMON.HEADER'
include 'COMMON.SPLITELE'
include 'COMMON.SHIELD'
include 'COMMON.GEO'
+ integer i
+ integer KDIAG,ICORFL,IXDR
COMMON /MACHSW/ KDIAG,ICORFL,IXDR
character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
character*80 ucase
character*320 controlcard
+ double precision seed
nglob_csa=0
eglob_csa=1d99
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
+ out_cart=index(controlcard,'OUT_CART').gt.0
+ out_int=index(controlcard,'OUT_INT').gt.0
+ gmatout=index(controlcard,'GMATOUT').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
out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
call readi(controlcard,'SYM',symetr,1)
+ call readi(controlcard,'PERMUT',npermut,1)
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)
- 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
+c call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
+c call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
+c call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
+c call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
+c call reada(controlcard,'DRMS',drms,0.1D0)
+c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+c write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
+c write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
+c write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
+c write (iout,'(a,f10.1)')'DRMS = ',drms
+cc write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
+c write (iout,'(a,f10.1)') 'Time limit (min):',timlim
+c endif
call readi(controlcard,'NZ_START',nz_start,0)
call readi(controlcard,'NZ_END',nz_end,0)
+ call readi(controlcard,'NRAN_START',nran_start,0)
+ write (iout,*) "nran_start",nran_start
c call readi(controlcard,'IZ_SC',iz_sc,0)
timlim=60.0D0*timlim
safety = 60.0d0*safety
- timem=timlim
modecalc=0
+ call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
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,'NOSEARCHSC').gt.0)
- searchsc=.not.searchsc
+ searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
+ & index(controlcard,'NOSEARCHSC').eq.0)
sideadd=(index(controlcard,'SIDEADD').gt.0)
energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
mremd_dec=(index(controlcard,'MREMD_DEC').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)
AFMlog=(index(controlcard,'AFM'))
selfguide=(index(controlcard,'SELFGUIDE'))
c print *,'AFMlog',AFMlog,selfguide,"KUPA"
indphi=index(controlcard,'PHI')
indback=index(controlcard,'BACK')
iranconf=index(controlcard,'RAND_CONF')
+ start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
+ extconf=(index(controlcard,'EXTCONF').gt.0)
+ if (start_from_model) then
+ iranconf=0
+ extconf=.false.
+ endif
i2ndstr=index(controlcard,'USE_SEC_PRED')
gradout=index(controlcard,'GRADOUT').gt.0
gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
call reada(controlcard,'BOXX',boxxsize,100.0d0)
call reada(controlcard,'BOXY',boxysize,100.0d0)
call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+ write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
c Cutoff range for interactions
- call reada(controlcard,"R_CUT",r_cut,15.0d0)
+ call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
+ call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+ write (iout,*) "Cutoff on interactions",r_cut_int
+ write (iout,*)
+ & "Cutoff in switching short and long range interactions in RESPA",
+ & r_cut_respa
+ write (iout,*) "lambda in switch function",rlamb
call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
if (lipthick.gt.0.0d0) then
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
C
C Read REMD settings
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.TIME1'
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.INTERACT'
include 'COMMON.NAMES'
include 'COMMON.GEO'
character*80 ucase
character*320 controlcard
character*3200 controlcard1
- integer iremd_m_total
+ integer iremd_m_total,i
if(me.eq.king.or..not.out1file)
& write (iout,*) "REMD setup"
C
C Read MD settings
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.TIME1'
include 'COMMON.MD'
+ include 'COMMON.QRESTR'
#ifndef LANG0
include 'COMMON.LANGEVIN'
#else
+#ifdef FIVEDIAG
+ include 'COMMON.LANGEVIN.lang0.5diag'
+#else
include 'COMMON.LANGEVIN.lang0'
#endif
+#endif
include 'COMMON.INTERACT'
include 'COMMON.NAMES'
include 'COMMON.GEO'
include 'COMMON.FFIELD'
character*80 ucase
character*320 controlcard
+ integer i
+ double precision eta
call card_concat(controlcard)
call readi(controlcard,"NSTEP",n_timestep,1000000)
call readi(controlcard,"NTWE",ntwe,100)
call readi(controlcard,"NTWX",ntwx,1000)
+ call readi(controlcard,"REST_FREQ",irest_freq,1000)
call reada(controlcard,"DT",d_time,1.0d-1)
call reada(controlcard,"DVMAX",dvmax,2.0d1)
call reada(controlcard,"DAMAX",damax,1.0d1)
& "Initial time step of numerical integration:",d_time,
& " natural units"
write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
+ write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
+ & " A"
+ write(iout,'(a60,i5)')"Frequency of updating interaction list",
+ & imatupdate
+ write(iout,'(a60,i5)')"Restart writing frequency",irest_freq
if (RESPA) then
write (iout,'(2a,i4,a)')
& "A-MTS algorithm used; initial time step for fast-varying",
& " short-range forces split into",ntime_split," steps."
write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
- & r_cut," lambda",rlamb
+ & r_cut_respa," lambda",rlamb
endif
write (iout,'(2a,f10.5)')
& "Maximum acceleration threshold to reduce the time step",
C
C Read molecular data.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
- integer error_msg
+ integer error_msg,ierror,ierr,ierrcode
#endif
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.SBRIDGE'
include 'COMMON.HEADER'
include 'COMMON.CONTROL'
+ include 'COMMON.SAXS'
include 'COMMON.DBASE'
include 'COMMON.THREAD'
include 'COMMON.CONTACTS'
character*256 pdbfile
character*400 weightcard
character*80 weightcard_t,ucase
- dimension itype_pdb(maxres)
+ integer itype_pdb(maxres)
common /pizda/ itype_pdb
logical seq_comp,fail
double precision energia(0:n_ene)
double precision secprob(3,maxdih_constr)
+ double precision co
+ double precision phihel,phibet,sigmahel,sigmabet
+ integer iti,nsi,maxsi
integer ilen
external ilen
- integer tperm
+ integer iperm,tperm
+ integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
+ double precision sumv
C
C Read PDB structure if applicable
C
do i=1,nres
itype(i)=rescode(i,sequence(i),iscode)
enddo
-C Assign initial virtual bond lengths
-c do i=2,nres
-c vbld(i)=vbl
-c vbld_inv(i)=vblinv
-c enddo
-c if (itype(1).eq.ntyp1) then
-c vbld(2)=vbld(2)/2
-c vbld_inv(2)=vbld_inv(2)*2
-c endif
-c if (itype(nres).eq.ntyp1) then
-c vbld(nres)=vbld(nres)/2
-c vbld_inv(nres)=vbld_inv(nres)*2
-c endif
-c do i=2,nres-1
-c vbld(i+nres)=dsc(iabs(itype(i)))
-c 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)
-c enddo
endif
c print *,nres
c print '(20i4)',(itype(i),i=1,nres)
do i=1,nres-1
write (iout,*) i,itype(i),itel(i)
enddo
- print *,'Call Read_Bridge.'
+c print *,'Call Read_Bridge.'
endif
nnt=1
nct=nres
cd print *,'NNT=',NNT,' NCT=',NCT
call seq2chains(nres,itype,nchain,chain_length,chain_border,
& ireschain)
+ chain_border1(1,1)=1
+ chain_border1(2,1)=chain_border(2,1)+1
+ do i=2,nchain-1
+ chain_border1(1,i)=chain_border(1,i)-1
+ chain_border1(2,i)=chain_border(2,i)+1
+ enddo
+ if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
+ chain_border1(2,nchain)=nres
write(iout,*) "nres",nres," nchain",nchain
do i=1,nchain
write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
- & chain_border(2,i)
+ & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
enddo
call chain_symmetry(nchain,nres,itype,chain_border,
- & chain_length,npermchain,tabpermchain)
+ & chain_length,npermchain,tabpermchain,nchain_group,nequiv,
+ & iequiv,mapchain)
c do i=1,nres
c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
c & ii=1,npermchain)
c enddo
+#ifdef DEBUG
write(iout,*) "residue permutations"
do i=1,nres
write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
enddo
+ call flush(iout)
+#endif
if (itype(1).eq.ntyp1) nnt=2
if (itype(nres).eq.ntyp1) nct=nct-1
+ write (iout,*) "nnt",nnt," nct",nct
+ call flush(iout)
#ifdef DFA
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!'
+c print*, 'init_dfa_vars finished!'
call read_dfa_info
- print*, 'read_dfa_info finished!'
+c print*, 'read_dfa_info finished!'
endif
#endif
if (pdbref) then
call reada(weightcard,"WDIHC",wdihc,0.591D0)
write (iout,*) "Weight of dihedral angle restraints",wdihc
read(inp,'(9x,3f7.3)')
+c & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
& (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
+ iconstr_dih=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
+ iconstr_dih(i)=ndih_constr
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
+ if (sumv.gt.0.0d0) then
+ do j=1,3
+ vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
+ enddo
+ else
+ vpsipred(1,ndih_constr)=1.0d0
+ vpsipred(2,ndih_constr)=0.0d0
+ vpsipred(3,ndih_constr)=0.0d0
+ endif
phibound(1,ndih_constr)=phihel*deg2rad
phibound(2,ndih_constr)=phibet*deg2rad
sdihed(1,ndih_constr)=sigmahel*deg2rad
endif
endif
c print *, "A TU"
- write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+c 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
+c write (iout,*) "After read_dist_constr nhpb",nhpb
if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
call hpb_partition
call NMRpeak_partition
endif
write (iout,*) "calling read_saxs_consrtr",nsaxs
if (nsaxs.gt.0) call read_saxs_constr
-
+c write (iout,*) "After read_saxs_constr"
+c call flush(iout)
if (constr_homology.gt.0) then
+c write (iout,*) "Calling read_constr_homology"
+c call flush(iout)
call read_constr_homology
if (indpdb.gt.0 .or. pdbref) then
do i=1,2*nres
enddo
else
homol_nset=0
+ if (start_from_model) then
+ nmodel_start=0
+ do
+ read(inp,'(a)',end=332,err=332) pdbfile
+ if (me.eq.king .or. .not. out1file)
+ & write (iout,'(a,5x,a)') 'Opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=336)
+ goto 335
+ 336 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ call flush(iout)
+ stop
+ 335 continue
+ unres_pdb=.false.
+ nres_temp=nres
+c call readpdb
+ call readpdb_template(nmodel_start+1)
+ close(ipdbin)
+ if (nres.ge.nres_temp) then
+ nmodel_start=nmodel_start+1
+ pdbfiles_chomo(nmodel_start)=pdbfile
+ do i=1,2*nres
+ do j=1,3
+ chomo(j,i,nmodel_start)=c(j,i)
+ enddo
+ enddo
+ else
+ if (me.eq.king .or. .not. out1file)
+ & write (iout,'(a,2i7,1x,a)')
+ & "Different number of residues",nres_temp,nres,
+ & " model skipped."
+ endif
+ nres=nres_temp
+ enddo
+ 332 continue
+ if (nmodel_start.eq.0) then
+ if (me.eq.king .or. .not. out1file)
+ & write (iout,'(a)')
+ & "No valid starting model found START_FROM_MODELS is OFF"
+ start_from_model=.false.
+ endif
+ write (iout,*) "nmodel_start",nmodel_start
+ endif
endif
& 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. .not. extconf .and. .not.
+ & start_from_model) 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)
- write (iout,*) "Exit READ_CART"
+ if (nnt.gt.1) c(:,nres+1)=c(:,1)
+ if (nct.lt.nres) c(:,2*nres)=c(:,nres)
+c write (iout,*) "Exit READ_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 j=1,3
+ dc(j,0)=c(j,1)
+ enddo
do i=1,nres-1
do j=1,3
dc(j,i)=c(j,i+1)-c(j,i)
enddo
call bond_regular
call chainbuild_extconf
- else
+ else if (.not. start_from_model) then
if(me.eq.king.or..not.out1file)
& write (iout,'(a)') 'Random-generated initial geometry.'
call bond_regular
c-----------------------------------------------------------------------------
subroutine read_bridge
C Read information about disulfide bridges.
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
+ integer ierror
#endif
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.THREAD'
include 'COMMON.TIME1'
include 'COMMON.SETUP'
+ integer i,j
C Read bridging residues.
read (inp,*) ns,(iss(i),i=1,ns)
- print *,'ns=',ns
+c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
+c numbers
+ icys=0
+ do i=1,ns
+ icys(iss(i))=i
+ enddo
+c 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.
enddo
write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
20 continue
- dhpb(i)=dbr
- forcon(i)=fbr
+c dhpb(i)=dbr
+c forcon(i)=fbr
enddo
do i=1,nss
ihpb(i)=ihpb(i)+nres
end
c----------------------------------------------------------------------------
subroutine read_x(kanal,*)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.CONTROL'
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
+ integer i,j,k,l,kanal
c Read coordinates from input
c
read(kanal,'(8f10.5)',end=10,err=10)
end
c----------------------------------------------------------------------------
subroutine read_threadbase
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.DBASE'
include 'COMMON.THREAD'
include 'COMMON.TIME1'
+ integer i,j,k
+ double precision dist
C Read pattern database for threading.
read (icbase,*) nseq
do i=1,nseq
end
c------------------------------------------------------------------------------
subroutine setup_var
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.DBASE'
include 'COMMON.THREAD'
include 'COMMON.TIME1'
+ integer i
C Set up variable list.
ntheta=nres-2
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.ntyp1) then
- nside=nside+1
+ nside=nside+1
ialph(i,1)=nvar+nside
- ialph(nside,2)=i
+ ialph(nside,2)=i
endif
enddo
if (indphi.gt.0) then
else
nvar=nvar+2*nside
endif
- 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)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.CONTROL'
include 'COMMON.DBASE'
include 'COMMON.THREAD'
+ include 'COMMON.SPLITELE'
include 'COMMON.TIME1'
- dimension itype_pdb(maxres)
+ integer i,j,itype_pdb(maxres)
common /pizda/ itype_pdb
+ double precision dd
+ double precision dist
character*2 iden
cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
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
+c 5/24/2020 Adam: Cutoff included to reduce array size
+ dd = dist(i,j)
+ if (dd.gt.r_cut_int) cycle
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)
+ dhpb(nhpb)=dd
enddo
enddo
cd write (iout,'(a)') 'Distance constraints:'
end
c----------------------------------------------------------------------------
subroutine map_read
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.MAP'
include 'COMMON.IOUNITS'
+ integer imap
character*3 angid(4) /'THE','PHI','ALP','OME'/
character*80 mapcard,ucase
do imap=1,nmap
end
c----------------------------------------------------------------------------
subroutine csaread
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
return
end
c----------------------------------------------------------------------------
-cfmc subroutine mcmfread
-cfmc implicit real*8 (a-h,o-z)
-cfmc include 'DIMENSIONS'
-cfmc include 'COMMON.MCMF'
-cfmc include 'COMMON.IOUNITS'
-cfmc include 'COMMON.GEO'
-cfmc character*80 ucase
-cfmc character*620 mcmcard
-cfmc call card_concat(mcmcard)
-cfmc
-cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
-cfmc write(iout,*)'MAXRANT=',maxrant
-cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
-cfmc write(iout,*)'MAXFAM=',maxfam
-cfmc call readi(mcmcard,'NNET1',nnet1,5)
-cfmc write(iout,*)'NNET1=',nnet1
-cfmc call readi(mcmcard,'NNET2',nnet2,4)
-cfmc write(iout,*)'NNET2=',nnet2
-cfmc call readi(mcmcard,'NNET3',nnet3,4)
-cfmc write(iout,*)'NNET3=',nnet3
-cfmc call readi(mcmcard,'ILASTT',ilastt,0)
-cfmc write(iout,*)'ILASTT=',ilastt
-cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
-cfmc write(iout,*)'MAXSTR=',maxstr
-cfmc maxstr_f=maxstr/maxfam
-cfmc write(iout,*)'MAXSTR_F=',maxstr_f
-cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
-cfmc write(iout,*)'NMCMF=',nmcmf
-cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
-cfmc write(iout,*)'IFOCUS=',ifocus
-cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
-cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
-cfmc call readi(mcmcard,'INTPRT',intprt,1000)
-cfmc write(iout,*)'INTPRT=',intprt
-cfmc call readi(mcmcard,'IPRT',iprt,100)
-cfmc write(iout,*)'IPRT=',iprt
-cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
-cfmc write(iout,*)'IMAXTR=',imaxtr
-cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
-cfmc write(iout,*)'MAXEVEN=',maxeven
-cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
-cfmc write(iout,*)'MAXEVEN1=',maxeven1
-cfmc call readi(mcmcard,'INIMIN',inimin,200)
-cfmc write(iout,*)'INIMIN=',inimin
-cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
-cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
-cfmc call readi(mcmcard,'NTHREAD',nthread,5)
-cfmc write(iout,*)'NTHREAD=',nthread
-cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
-cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
-cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
-cfmc write(iout,*)'MAXPERT=',maxpert
-cfmc call readi(mcmcard,'IRMSD',irmsd,1)
-cfmc write(iout,*)'IRMSD=',irmsd
-cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
-cfmc write(iout,*)'DENEMIN=',denemin
-cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
-cfmc write(iout,*)'RCUT1S=',rcut1s
-cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
-cfmc write(iout,*)'RCUT1E=',rcut1e
-cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
-cfmc write(iout,*)'RCUT2S=',rcut2s
-cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
-cfmc write(iout,*)'RCUT2E=',rcut2e
-cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
-cfmc write(iout,*)'DPERT1=',d_pert1
-cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
-cfmc write(iout,*)'DPERT1A=',d_pert1a
-cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
-cfmc write(iout,*)'DPERT2=',d_pert2
-cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
-cfmc write(iout,*)'DPERT2A=',d_pert2a
-cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
-cfmc write(iout,*)'DPERT2B=',d_pert2b
-cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
-cfmc write(iout,*)'DPERT2C=',d_pert2c
-cfmc d_pert1=deg2rad*d_pert1
-cfmc d_pert1a=deg2rad*d_pert1a
-cfmc d_pert2=deg2rad*d_pert2
-cfmc d_pert2a=deg2rad*d_pert2a
-cfmc d_pert2b=deg2rad*d_pert2b
-cfmc d_pert2c=deg2rad*d_pert2c
-cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
-cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
-cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
-cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
-cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
-cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
-cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
-cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
-cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
-cfmc write(iout,*)'RCUTINI=',rcutini
-cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
-cfmc write(iout,*)'GRAT=',grat
-cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
-cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
-cfmc
-cfmc return
-cfmc end
-c----------------------------------------------------------------------------
subroutine mcmread
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.MCM'
include 'COMMON.MCE'
include 'COMMON.IOUNITS'
character*80 ucase
character*320 mcmcard
+ integer i
call card_concat(mcmcard)
call readi(mcmcard,'MAXACC',maxacc,100)
call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
end
c----------------------------------------------------------------------------
subroutine read_minim
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.MINIM'
include 'COMMON.IOUNITS'
end
c----------------------------------------------------------------------------
subroutine read_angles(kanal,*)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
+ integer i,kanal
c Read angles from input
c
read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
end
c----------------------------------------------------------------------------
subroutine openunits
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
+ integer ierror
character*16 form,nodename
integer nodelen
#endif
include 'COMMON.IOUNITS'
include 'COMMON.MD'
include 'COMMON.CONTROL'
- integer lenpre,lenpot,ilen,lentmp
+ integer lenpre,lenpot,ilen,lentmp,npos
external ilen
character*3 out1file_text,ucase
character*3 ll
#else
if (me.eq.king .or. .not.out1file)
& open(iout,file=outname,status='unknown')
+#define DEBUG
#ifdef DEBUG
if (fg_rank.gt.0) then
write (liczba,'(i3.3)') myrank/nfgtasks
& status='unknown')
endif
#endif
+#undef DEBUG
if(me.eq.king) then
open(igeom,file=intname,status='unknown',access='append')
open(ipdb,file=pdbname,status='unknown')
card=card(:ilen(card)+1)//karta
return
end
-c----------------------------------------------------------------------------------
+c------------------------------------------------------------------------------
subroutine readrst
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
include 'COMMON.MD'
+ include 'COMMON.QRESTR'
+#ifdef FIVEDIAG
+ include 'COMMON.LAGRANGE.5diag'
+#else
+ include 'COMMON.LAGRANGE'
+#endif
+ integer i,j
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
totTafm=totT
- do i=1,2*nres
+ do i=0,2*nres-1
read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
enddo
- do i=1,2*nres
+ do i=0,2*nres-1
read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
if(usampl) then
close(irest2)
return
end
-c---------------------------------------------------------------------------------
+c------------------------------------------------------------------------------
subroutine read_fragments
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.MD'
+ include 'COMMON.QRESTR'
include 'COMMON.CONTROL'
+ integer i
read(inp,*) nset,nfrag,npair,nfrag_back
loc_qlike=(nfrag_back.lt.0)
nfrag_back=iabs(nfrag_back)
end
C---------------------------------------------------------------------------
subroutine read_afminp
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.IOUNITS'
include 'COMMON.SBRIDGE'
character*320 afmcard
- print *, "wchodze"
+ integer i
+c 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
+c 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
+c print *,'initdist',distafminit
return
end
c-------------------------------------------------------------------------------
subroutine read_saxs_constr
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
#endif
include 'COMMON.SETUP'
include 'COMMON.CONTROL'
+ include 'COMMON.SAXS'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.SBRIDGE'
- double precision cm(3)
+ double precision cm(3),cnorm
+ integer i,j
c read(inp,*) nsaxs
write (iout,*) "Calling read_saxs nsaxs",nsaxs
call flush(iout)
c-------------------------------------------------------------------------------
subroutine read_dist_constr
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
include 'COMMON.IOUNITS'
include 'COMMON.SBRIDGE'
include 'COMMON.INTERACT'
- integer ifrag_(2,100),ipair_(2,1000)
+ integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
+ integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
double precision wfrag_(100),wpair_(1000)
+ double precision ddjk,dist,dist_cut,fordepthmax
character*5000 controlcard
logical normalize,next
integer restr_type
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+ & write (iout,'(a,3i6,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.restr ",
+ write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
enddo
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+ & write (iout,'(a,3i6,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.restr ",
+ write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
enddo
ipeak(2,npeak)=i
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & write (iout,'(a,5i6,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 ",
+ write (iout,'(a,5i6,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),
irestr_type(nhpb)=11
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & write (iout,'(a,4i6,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 ",
+ write (iout,'(a,4i6,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
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ & write (iout,'(a,4i6,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,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,4i6,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
#ifdef MPI
if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+ & write (iout,'(a,4i6,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 ",
+ write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
#endif
endif
dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
#ifdef MPI
if (.not.out1file .or. me.eq.king) then
- write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,4i6,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
- write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ write (iout,'(a,4i6,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)
end
c-------------------------------------------------------------------------------
subroutine read_constr_homology
-
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
#endif
include 'COMMON.SETUP'
include 'COMMON.CONTROL'
+ include 'COMMON.HOMOLOGY'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.MD'
+ include 'COMMON.QRESTR'
include 'COMMON.GEO'
include 'COMMON.INTERACT'
include 'COMMON.NAMES'
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
+ integer ki,i,ii,j,k,l,ii_in_use(maxdim_cont),i_tmp,idomain_tmp,
+ & irec,ik,iistart,nres_temp
integer ilen
external ilen
- logical liiflag
+ logical liiflag,lfirst
+ integer i01,i10
c
c FP - Nov. 2014 Temporary specifications for new vars
c
double precision, dimension (max_template,maxres) :: rescore
double precision, dimension (max_template,maxres) :: rescore2
double precision, dimension (max_template,maxres) :: rescore3
- character*24 pdbfile,tpl_k_rescore
+ double precision distal
+ character*24 tpl_k_rescore
+ character*256 pdbfile
c -----------------------------------------------------------------
c Reading multiple PDB ref structures and calculation of retraints
c not using pre-computed ones stored in files model_ki_{dist,angle}
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.
+ iranconf=(indpdb.le.0)
endif
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
- 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
+c if(start_from_model .and. rest) then
+c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+c write(iout,*) 'START_FROM_MODELS is OFF'
+c write(iout,*) 'remove restart keyword from input'
+c endif
+c endif
+ if (start_from_model) nmodel_start=constr_homology
if (homol_nset.gt.1)then
call card_concat(controlcard)
read(controlcard,*) (waga_homology(i),i=1,homol_nset)
if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
- write(iout,*) "iset homology_weight "
+c write(iout,*) "iset homology_weight "
do i=1,homol_nset
write(iout,*) i,waga_homology(i)
enddo
cd write (iout,*) "nnt",nnt," nct",nct
cd call flush(iout)
+ if (read_homol_frag) then
+ call read_klapaucjusz
+ else
lim_odl=0
lim_dih=0
c
c write(iout,*) 'nnt=',nnt,'nct=',nct
c
- do i = nnt,nct
- do k=1,constr_homology
- idomain(k,i)=0
- enddo
- enddo
-
- ii=0
- do i = nnt,nct-2
- do j=i+2,nct
- ii=ii+1
- ii_in_use(ii)=0
- enddo
- enddo
-
- if (read_homol_frag) then
- call read_klapaucjusz
- else
+c do i = nnt,nct
+c do k=1,constr_homology
+c idomain(k,i)=0
+c enddo
+c enddo
+ idomain=0
+c ii=0
+c do i = nnt,nct-2
+c do j=i+2,nct
+c ii=ii+1
+c ii_in_use(ii)=0
+c enddo
+c enddo
+ ii_in_use=0
do k=1,constr_homology
read(inp,'(a)') pdbfile
+ pdbfiles_chomo(k)=pdbfile
if(me.eq.king .or. .not. out1file)
& write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
& pdbfile(:ilen(pdbfile))
tpl_k_rescore="template"//kic2//".sco"
unres_pdb=.false.
+ nres_temp=nres
if (read2sigma) then
call readpdb_template(k)
else
call readpdb
endif
+ nres_chomo(k)=nres
+ nres=nres_temp
c
c Distance restraints
c
c ... --> odl(k,ii)
C Copy the coordinates from reference coordinates (?)
- do i=1,2*nres
+ do i=1,2*nres_chomo(k)
do j=1,3
c(j,i)=cref(j,i)
c write (iout,*) "c(",j,i,") =",c(j,i)
endif
sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
else
- ii=ii+1
- l_homo(k,ii)=.false.
+c ii=ii+1
+c l_homo(k,ii)=.false.
endif
enddo
enddo
lim_odl=ii
endif
+c write (iout,*) "Distance restraints set"
+c call flush(iout)
c
c Theta, dihedral and SC retraints
c
enddo
lim_dih=nct-nnt-2
endif
+c write (iout,*) "Dihedral angle restraints set"
+c call flush(iout)
if (waga_theta.gt.0.0d0) then
c open (ientin,file=tpl_k_sigma_theta,status='old')
c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
enddo
endif
+c write (iout,*) "Angle restraints set"
+c call flush(iout)
if (waga_d.gt.0.0d0) then
c open (ientin,file=tpl_k_sigma_d,status='old')
enddo
endif
enddo
+c write (iout,*) "SC restraints set"
+c call flush(iout)
c
c remove distance restraints not used in any model from the list
c shift data in all arrays
c
+c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
if (waga_dist.ne.0.0d0) then
ii=0
liiflag=.true.
+ lfirst=.true.
do i=nnt,nct-2
do j=i+2,nct
ii=ii+1
- if (ii_in_use(ii).eq.0.and.liiflag) then
+c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+c & .and. distal.le.dist2_cut ) then
+c write (iout,*) "i",i," j",j," ii",ii
+c call flush(iout)
+ if (ii_in_use(ii).eq.0.and.liiflag.or.
+ & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
liiflag=.false.
- iistart=ii
+ i10=ii
+ if (lfirst) then
+ lfirst=.false.
+ iistart=ii
+ else
+ if(i10.eq.lim_odl) i10=i10+1
+ do ki=0,i10-i01-1
+ ires_homo(iistart+ki)=ires_homo(ki+i01)
+ jres_homo(iistart+ki)=jres_homo(ki+i01)
+ ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+ do k=1,constr_homology
+ odl(k,iistart+ki)=odl(k,ki+i01)
+ sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+ l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+ enddo
+ enddo
+ iistart=iistart+i10-i01
+ endif
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
+ if (ii_in_use(ii).ne.0.and..not.liiflag) then
+ i01=ii
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
+ lim_odl=iistart-1
endif
-
+c write (iout,*) "Removing distances completed"
+c call flush(iout)
endif ! .not. klapaucjusz
if (constr_homology.gt.0) call homology_partition
+c write (iout,*) "After homology_partition"
+c call flush(iout)
if (constr_homology.gt.0) call init_int_table
-c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
-c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c write (iout,*) "After init_int_table"
+c call flush(iout)
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
c
c Print restraints
c
if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write (iout,*) "Distance restraints from templates"
do ii=1,lim_odl
- write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
+ write(iout,'(3i7,100(2f8.2,1x,l1,4x))')
& ii,ires_homo(ii),jres_homo(ii),
& (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
& ki=1,constr_homology)
enddo
write (iout,*) "Dihedral angle restraints from templates"
do i=nnt+3,nct
- write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+ write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
& (rad2deg*dih(ki,i),
& rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
enddo
write (iout,*) "Virtual-bond angle restraints from templates"
do i=nnt+2,nct
- write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+ write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)),
& (rad2deg*thetatpl(ki,i),
& rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
enddo
write (iout,*) "SC restraints from templates"
do i=nnt,nct
- write(iout,'(i5,100(4f8.2,4x))') i,
+ write(iout,'(i7,100(4f8.2,4x))') i,
& (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
& 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
enddo
#endif
c------------------------------------------------------------------------------
subroutine copy_to_tmp(source)
+ implicit none
include "DIMENSIONS"
include "COMMON.IOUNITS"
character*(*) source
end
c------------------------------------------------------------------------------
subroutine move_from_tmp(source)
+ implicit none
include "DIMENSIONS"
include "COMMON.IOUNITS"
character*(*) source
C
C Initialize random number generator
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
logical OKRandom, prng_restart
real*8 r1
integer iseed_array(4)
+ integer error_msg,ierr
#endif
include 'COMMON.IOUNITS'
include 'COMMON.TIME1'
include 'COMMON.MD'
include 'COMMON.FFIELD'
include 'COMMON.SETUP'
+ integer i,iseed
+ double precision seed,ran_number
iseed=-dint(dabs(seed))
if (iseed.eq.0) then
write (iout,'(/80(1h*)/20x,a/80(1h*))')
end
c----------------------------------------------------------------------
subroutine read_klapaucjusz
-
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include 'mpif.h'
#endif
include 'COMMON.SETUP'
include 'COMMON.CONTROL'
+ include 'COMMON.HOMOLOGY'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.MD'
include 'COMMON.NAMES'
character*256 fragfile
integer ninclust(maxclust),inclust(max_template,maxclust),
- & nresclust(maxclust),iresclust(maxres,maxclust)
+ & nresclust(maxclust),iresclust(maxres,maxclust),nclust
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
+ integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp,
+ & idomain_tmp,
+ & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,
+ & i01,i10,nnt_chain,nct_chain
+ integer itype_temp(maxres)
+ double precision distal
logical lprn /.true./
+ integer nres_temp
integer ilen
external ilen
- logical liiflag
+ logical liiflag,lfirst
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
+ character*24 tpl_k_rescore
+ character*256 pdbfile
c
c For new homol impl
c
include 'COMMON.VAR'
c
+c write (iout,*) "READ_KLAPAUCJUSZ"
+c print *,"READ_KLAPAUCJUSZ"
+c call flush(iout)
call getenv("FRAGFILE",fragfile)
+ write (iout,*) "Opening", fragfile
+ call flush(iout)
open(ientin,file=fragfile,status="old",err=10)
- read(ientin,*) constr_homology,nclust
- l_homo = .false.
+c write (iout,*) " opened"
+c call flush(iout)
+
sigma_theta=0.0
sigma_d=0.0
sigma_dih=0.0
+ l_homo = .false.
+
+ nres_temp=nres
+ itype_temp=itype
+ ii=0
+ lim_odl=0
+
+c write (iout,*) "Entering loop"
+c call flush(iout)
+
+ DO IGR = 1,NCHAIN_GROUP
+
+c write (iout,*) "igr",igr
+ call flush(iout)
+ read(ientin,*) constr_homology,nclust
+
+ if (start_from_model) then
+ nmodel_start=constr_homology
+ else
+ nmodel_start=0
+ endif
+
+ ii_old=lim_odl
+
+ ichain=iequiv(1,igr)
+ nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
+ nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
+c write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
+
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))
+ pdbfiles_chomo(k)=pdbfile
open(ipdbin,file=pdbfile,status='old',err=33)
goto 34
33 write (iout,'(a,5x,a)') 'Error opening PDB file',
stop
34 continue
unres_pdb=.false.
+c nres_temp=nres
call readpdb_template(k)
+ nres_chomo(k)=nres
+c nres=nres_temp
do i=1,nres
rescore(k,i)=0.2d0
rescore2(k,i)=1.0d0
do ll = 1,ninclust(l)
k = inclust(ll,l)
+c write (iout,*) "l",l," ll",ll," k",k
do i=1,nres
idomain(k,i)=0
enddo
c
c ... --> odl(k,ii)
C Copy the coordinates from reference coordinates (?)
+c nres_temp=nres
+ nres=nres_chomo(k)
do i=1,2*nres
do j=1,3
c(j,i)=chomo(j,i,k)
enddo
enddo
call int_from_cart(.true.,.false.)
- call sc_loc_geom(.true.)
+ call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
enddo
+c nres=nres_temp
if (waga_dist.ne.0.0d0) then
- ii=0
- do i = nnt,nct-2
- do j=i+2,nct
+ ii=ii_old
+c do i = nnt,nct-2
+ do i = nnt_chain,nct_chain-2
+c do j=i+2,nct
+ do j=i+2,nct_chain
x12=c(1,i)-c(1,j)
y12=c(2,i)-c(2,j)
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
+c & constr_homology
+ ires_homo(ii)=i+chain_border1(1,igr)-1
+ jres_homo(ii)=j+chain_border1(1,igr)-1
odl(k,ii)=distal
if (read2sigma) then
sigma_odl(k,ii)=0
c Theta, dihedral and SC retraints
c
if (waga_angle.gt.0.0d0) then
- do i = nnt+3,nct
+ do i = nnt_chain+3,nct_chain
+ iii=i+chain_border1(1,igr)-1
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)+
+ dih(k,iii)=phiref(i)
+ sigma_dih(k,iii)=
+ & (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))
+ if (sigma_dih(k,iii).ne.0)
+ & sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
enddo
- lim_dih=nct-nnt-2
+c lim_dih=nct-nnt-2
endif
if (waga_theta.gt.0.0d0) then
- do i = nnt+2,nct
+ do i = nnt_chain+2,nct_chain
+ iii=i+chain_border1(1,igr)-1
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)+
+ thetatpl(k,iii)=thetaref(i)
+ sigma_theta(k,iii)=(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))
+ if (sigma_theta(k,iii).ne.0)
+ & sigma_theta(k,iii)=1.0d0/
+ & (sigma_theta(k,iii)*sigma_theta(k,iii))
enddo
endif
if (waga_d.gt.0.0d0) then
- do i = nnt,nct
+ do i = nnt_chain,nct_chain
+ iii=i+chain_border1(1,igr)-1
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
+ xxtpl(k,iii)=xxref(i)
+ yytpl(k,iii)=yyref(i)
+ zztpl(k,iii)=zzref(i)
+ sigma_d(k,iii)=rescore(k,i)
+ if (sigma_d(k,iii).ne.0)
+ & sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
+c if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
enddo
endif
enddo ! l
c remove distance restraints not used in any model from the list
c shift data in all arrays
c
+c write (iout,*) "ii_old",ii_old
if (waga_dist.ne.0.0d0) then
- ii=0
+#ifdef DEBUG
+ write (iout,*) "Distance restraints from templates"
+ do iii=1,lim_odl
+ write(iout,'(4i5,100(2f8.2,1x,l1,4x))')
+ & iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii),
+ & (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii),
+ & ki=1,constr_homology)
+ enddo
+#endif
+ ii=ii_old
liiflag=.true.
- do i=nnt,nct-2
- do j=i+2,nct
+ lfirst=.true.
+ do i=nnt_chain,nct_chain-2
+ do j=i+2,nct_chain
ii=ii+1
- if (ii_in_use(ii).eq.0.and.liiflag) then
+c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+c & .and. distal.le.dist2_cut ) then
+c write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
+c call flush(iout)
+ if (ii_in_use(ii).eq.0.and.liiflag.or.
+ & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
liiflag=.false.
- iistart=ii
+ i10=ii
+ if (lfirst) then
+ lfirst=.false.
+ iistart=ii
+ else
+ if(i10.eq.lim_odl) i10=i10+1
+ do ki=0,i10-i01-1
+ ires_homo(iistart+ki)=ires_homo(ki+i01)
+ jres_homo(iistart+ki)=jres_homo(ki+i01)
+ ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+ do k=1,constr_homology
+ odl(k,iistart+ki)=odl(k,ki+i01)
+ sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+ l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+ enddo
+ enddo
+ iistart=iistart+i10-i01
+ endif
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
+ if (ii_in_use(ii).ne.0.and..not.liiflag) then
+ i01=ii
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
+ lim_odl=iistart-1
endif
+ lll=lim_odl-ii_old
+
+ do i=2,nequiv(igr)
+
+ ichain=iequiv(i,igr)
+
+ do j=nnt_chain,nct_chain
+ jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
+ do k=1,constr_homology
+ dih(k,jj)=dih(k,j)
+ sigma_dih(k,jj)=sigma_dih(k,j)
+ thetatpl(k,jj)=thetatpl(k,j)
+ sigma_theta(k,jj)=sigma_theta(k,j)
+ xxtpl(k,jj)=xxtpl(k,j)
+ yytpl(k,jj)=yytpl(k,j)
+ zztpl(k,jj)=zztpl(k,j)
+ sigma_d(k,jj)=sigma_d(k,j)
+ enddo
+ enddo
+
+ jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
+c write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
+ do j=ii_old+1,lim_odl
+ ires_homo(j+lll)=ires_homo(j)+jj
+ jres_homo(j+lll)=jres_homo(j)+jj
+ do k=1,constr_homology
+ odl(k,j+lll)=odl(k,j)
+ sigma_odl(k,j+lll)=sigma_odl(k,j)
+ l_homo(k,j+lll)=l_homo(k,j)
+ enddo
+ enddo
+
+ ii_old=ii_old+lll
+ lim_odl=lim_odl+lll
+
+ enddo
+
+ ENDDO ! IGR
+
+ if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
+ nres=nres_temp
+ itype=itype_temp
+
return
10 stop "Error in fragment file"
end
projects / unres.git / blobdiff