X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Freadrtns_CSA.F;h=d76b29e0bf84b2232c41dc8c0975b47ccd4ce87a;hb=1441ec2b7d68b4fd94014649e4bd9574db0bfc36;hp=66f7c17924c5364ba6ff6ab2f29b9a981a43436a;hpb=020e579626d686ec20ecd9f0cc4c8313f474e152;p=unres.git diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index 66f7c17..d76b29e 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -1,5 +1,5 @@ subroutine readrtns - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -9,6 +9,7 @@ include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' include 'COMMON.SPLITELE' + integer i,j logical file_exist C Read job setup parameters call read_control @@ -84,18 +85,19 @@ C------------------------------------------------------------------------------- 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' @@ -109,10 +111,13 @@ C 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 @@ -125,6 +130,9 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file 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 @@ -144,33 +152,33 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword 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) 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) @@ -179,7 +187,6 @@ c call readi(controlcard,'IZ_SC',iz_sc,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" @@ -288,6 +295,12 @@ cfmc modecalc=10 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 @@ -297,9 +310,16 @@ 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) + 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 @@ -328,25 +348,6 @@ C endif 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 @@ -360,7 +361,7 @@ c-------------------------------------------------------------------------- C C Read REMD settings C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.TIME1' @@ -368,8 +369,12 @@ C #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' @@ -379,7 +384,7 @@ C 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" @@ -445,16 +450,21 @@ c-------------------------------------------------------------------------- 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' @@ -464,6 +474,8 @@ C 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) @@ -540,12 +552,16 @@ c if performing umbrella sampling, fragments constrained are read from the frag & "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 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", @@ -681,11 +697,11 @@ c------------------------------------------------------------------------------ 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' @@ -698,6 +714,7 @@ C include 'COMMON.SBRIDGE' include 'COMMON.HEADER' include 'COMMON.CONTROL' + include 'COMMON.SAXS' include 'COMMON.DBASE' include 'COMMON.THREAD' include 'COMMON.CONTACTS' @@ -713,14 +730,19 @@ C 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 @@ -789,25 +811,6 @@ C Convert sequence to numeric code 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) @@ -833,17 +836,25 @@ 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) @@ -855,15 +866,18 @@ c enddo do i=1,nres write(iout,*) i,(iperm(i,ii),ii=1,npermchain) enddo + call flush(iout) 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 @@ -1088,10 +1102,10 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) 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 @@ -1154,6 +1168,49 @@ c print *, "A TU" 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 + call readpdb + 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,2i5,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 @@ -1163,14 +1220,15 @@ C 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" +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) @@ -1372,10 +1430,11 @@ c-------------------------------------------------------------------------- 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' @@ -1392,9 +1451,16 @@ C Read information about disulfide bridges. 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. @@ -1447,8 +1513,8 @@ C bridging residues. 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 @@ -1460,7 +1526,7 @@ C bridging residues. end c---------------------------------------------------------------------------- subroutine read_x(kanal,*) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -1469,6 +1535,7 @@ c---------------------------------------------------------------------------- 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) @@ -1499,7 +1566,7 @@ c end c---------------------------------------------------------------------------- subroutine read_threadbase - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' @@ -1515,6 +1582,8 @@ c---------------------------------------------------------------------------- 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 @@ -1544,7 +1613,7 @@ c & nres_base(1,i)) end c------------------------------------------------------------------------------ subroutine setup_var - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' @@ -1560,17 +1629,17 @@ c------------------------------------------------------------------------------ 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 @@ -1580,13 +1649,12 @@ C Set up variable list. 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' @@ -1601,9 +1669,12 @@ C Generate CA distance constraints. 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, @@ -1613,11 +1684,14 @@ cd & ' nsup',nsup 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:' @@ -1638,10 +1712,11 @@ cd enddo 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 @@ -1684,7 +1759,7 @@ c---------------------------------------------------------------------------- end c---------------------------------------------------------------------------- subroutine csaread - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.GEO' @@ -1746,114 +1821,15 @@ c!bankt 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) @@ -1913,7 +1889,7 @@ C Probabilities of different move types end c---------------------------------------------------------------------------- subroutine read_minim - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.MINIM' include 'COMMON.IOUNITS' @@ -1939,13 +1915,14 @@ c---------------------------------------------------------------------------- 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) @@ -2039,10 +2016,11 @@ c---------------------------------------------------------------------------- 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 @@ -2050,7 +2028,7 @@ c---------------------------------------------------------------------------- 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 @@ -2439,20 +2417,28 @@ c---------------------------------------------------------------------------- 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 @@ -2461,9 +2447,9 @@ c------------------------------------------------------------------------------- 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' @@ -2472,7 +2458,9 @@ c------------------------------------------------------------------------------- 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) @@ -2522,7 +2510,7 @@ c------------------------------------------------------------------------------- end C--------------------------------------------------------------------------- subroutine read_afminp - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -2533,35 +2521,38 @@ C--------------------------------------------------------------------------- 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) @@ -2618,7 +2609,7 @@ c SAXS "spheres". c------------------------------------------------------------------------------- subroutine read_dist_constr - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -2629,8 +2620,10 @@ c------------------------------------------------------------------------------- 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 @@ -3017,16 +3010,18 @@ C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(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' @@ -3043,10 +3038,12 @@ c & sigma_odl_temp(maxres,maxres,max_template) 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),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 @@ -3055,6 +3052,7 @@ c double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 double precision, dimension (max_template,maxres) :: rescore3 + double precision distal character*24 pdbfile,tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -3078,15 +3076,17 @@ c Alternative: reading from input 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) @@ -3151,17 +3151,20 @@ c 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) @@ -3251,6 +3254,8 @@ c & constr_homology enddo lim_odl=ii endif +c write (iout,*) "Distance restraints set" +c call flush(iout) c c Theta, dihedral and SC retraints c @@ -3291,6 +3296,8 @@ c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) 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') @@ -3326,6 +3333,8 @@ c rescore(k,i-2) ! right expression ? 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') @@ -3361,51 +3370,66 @@ c read (ientin,*) sigma_d(k,i) ! 1st variant 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 @@ -3455,6 +3479,7 @@ c---------------------------------------------------------------------- #endif c------------------------------------------------------------------------------ subroutine copy_to_tmp(source) + implicit none include "DIMENSIONS" include "COMMON.IOUNITS" character*(*) source @@ -3474,6 +3499,7 @@ c------------------------------------------------------------------------------ end c------------------------------------------------------------------------------ subroutine move_from_tmp(source) + implicit none include "DIMENSIONS" include "COMMON.IOUNITS" character*(*) source @@ -3490,13 +3516,14 @@ c------------------------------------------------------------------------------ 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' @@ -3512,6 +3539,8 @@ C 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*))') @@ -3573,13 +3602,14 @@ c r1 = prng_next(me) 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' @@ -3588,13 +3618,16 @@ c---------------------------------------------------------------------- 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, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp, + & ik,ll,ii,kk,iistart,iishift,lim_xx + double precision distal logical lprn /.true./ + integer nres_temp integer ilen external ilen logical liiflag @@ -3629,7 +3662,10 @@ c Read pdb files stop 34 continue unres_pdb=.false. + nres_temp=nres call readpdb_template(k) + nres_chomo(k)=nres + nres=nres_temp do i=1,nres rescore(k,i)=0.2d0 rescore2(k,i)=1.0d0 @@ -3663,6 +3699,8 @@ c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) + nres_temp=nres + nres=nres_chomo(k) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) @@ -3670,11 +3708,12 @@ c write (iout,*) "c(",j,i,") =",c(j,i) 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 + nres=nres_temp if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2