subroutine read_general_data(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,kolor,klucz #endif include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.VMCPAR" include "COMMON.TORSION" include "COMMON.INTERACT" include "COMMON.ENERGIES" include "COMMON.MINPAR" include "COMMON.IOUNITS" include "COMMON.TIME1" include "COMMON.PROTFILES" include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.THERMAL" include "COMMON.FFIELD" include "COMMON.OPTIM" include "COMMON.CONTROL" include "COMMON.SCCOR" include "COMMON.SPLITELE" character*800 controlcard integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le integer ind,itype1,itype2,itypf,itypsc,itypp,itypt1,itypt2 integer ilen,lenpot,lenpre external ilen character*4 liczba,liczba1 #ifndef ISNAN external proc_proc #endif character*16 ucase character*16 key external ucase double precision xchuj,weitemp,weitemp_low,weitemp_up character*80 item(7) integer nitem integer rescode c write (iout,*) "Enter read_general_data" c call flush(iout) lenpot=ilen(pot) lenpre=ilen(prefix) do i=1,max_ene ename(i)=" " enddo C Read first record: seed and number of energy components call card_concat(controlcard,.true.) c write (iout,*) "card1",controlcard c call flush(iout) call readi(controlcard,"ISEED",iseed,0) if (iseed.eq.0) stop "Seed is zero!" c print *,me," iseed",iseed call readi(controlcard,"NPARMSET",nparmset,1) c print *,me," nparmset",nparmset #ifdef MPI c Split processor pool if multiple parameter sets are treated if (nparmset.eq.1) then WHAM_COMM = MPI_COMM_WORLD c print *,me," opening ",outname," opened" open(iout,file=outname,status='unknown') myparm=1 c print *,me," outname ",outname," opened" else if (nprocs.lt.nparmset) then write (*,*) & "*** Cannot split parameter sets for fewer processors than sets", & nprocs,nparmset call MPI_Finalize(ierror) stop endif c write (iout,*) "nparmset",nparmset nprocs = nprocs/nparmset kolor = me/nprocs klucz = mod(me,nprocs) c write (*,*) "My old rank",me," kolor",kolor," klucz",klucz call MPI_Comm_split(MPI_COMM_WORLD,kolor,klucz,WHAM_COMM,ierror) call MPI_Comm_size(WHAM_COMM,nprocs,ierror) call MPI_Comm_rank(WHAM_COMM,me,ierror) c write (*,*) "My new rank",me," comm size",nprocs c write (*,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, c & " WHAM_COMM",WHAM_COMM myparm=kolor+1 c write (*,*) "My parameter set is",myparm write(liczba,'(bz,i4.4)') me write(liczba1,'(bz,i4.4)') myparm outname=prefix(:lenpre)//'.out_par'//liczba1//'_'// & pot(:lenpot)//liczba open(iout,file=outname,status='unknown') endif #endif c print *,me," checkpoint 1" energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) call readi(controlcard,'TORMODE',tor_mode,0) c print *,me," tor_mode",tor_mode call readi(controlcard,'SHIELD',shield_mode,0) call readi(controlcard,"N_ENE",n_ene,max_ene) c print *,"iseed",iseed," n_ene",n_ene call readi(controlcard,"NPARMSET",nparmset,1) geom_and_wham_weights = & index(controlcard,"GEOM_AND_WHAM_WEIGHTS").gt.0 c write (iout,*) "GEOM_AND_WHAM_WEIGHTS",geom_and_wham_weights if (iseed.gt.0) iseed=-iseed call vrndst(iseed) out_newe=index(controlcard,"OUT_NEWE").gt.0 unres_pdb = index(controlcard,"UNRES_PDB").gt.0 c write (iout,*) "UNRES_PDB ",unres_pdb c Energy calculation settings section call readi(controlcard,'TORMODE',tor_mode,0) C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) c print *,me," checkpoint 2" c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) call reada(controlcard,"LIPTHICK",lipthick,0.0d0) call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) if (lipthick.gt.0.0d0) then bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick C endif if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick if ((lipbufthick*2.0d0).gt.lipthick) &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" endif c write(iout,*) "bordliptop=",bordliptop c write(iout,*) "bordlipbot=",bordlipbot c write(iout,*) "bufliptop=",bufliptop c write(iout,*) "buflipbot=",buflipbot call readi(controlcard,'SYM',symetr,1) c print *,me," checkpoint 3" c write (iout,*) "n_ene",n_ene call flush(iout) c Energy-term-weights section wname(4)="WCORRH" C Read third record: weights do i=1,max_ene ww(i)=ww0(i) enddo c print *,me," checkpoint 4" C Read fourth record: masks call card_concat(controlcard,.true.) c write (iout,*) "card2",controlcard do i=1,n_ene key = "MASK_"//wname(i)(2:ilen(wname(i))) call readi(controlcard,key(:ilen(key)),imask(i),0) enddo c print *,me," checkpoint 5" C Read fifth record: lower limits of weights call card_concat(controlcard,.true.) c write (iout,*) "card3",controlcard do i=1,n_ene key = "WLOW_"//wname(i)(2:ilen(wname(i))) call reada(controlcard,key(:ilen(key)),ww_low(i),ww(i)) enddo C Read sixth record: upper limits of weights call card_concat(controlcard,.true.) c write (iout,*) "card4",controlcard do i=1,n_ene key = "WUP_"//wname(i)(2:ilen(wname(i))) call reada(controlcard,key(:ilen(key)),ww_up(i),ww(i)) enddo C Read seventh record: VMC parameters call card_concat(controlcard,.true.) c write (iout,*) "card5",controlcard call readi(controlcard,"MODE",mode,3) call readi(controlcard,"NSCANCYCLE",nscancycle,3) call readi(controlcard,"MAXSTEP_SCAN",maxstep_scan,50) call reada(controlcard,"RSTEP_SCAN",step_scan,1.0d-1) call readi(controlcard,"READ_STAT",read_stat,3) call readi(controlcard,"RESCALE_MODE",rescale_mode,2) c init_ene = index(controlcard,"INIT_ENE").gt.0 .and. read_stat.gt.1 init_ene = .true. call readi(controlcard,"NMCM",nmcm,0) call readi(controlcard,"MAXSCAN",maxscan,0) call readi(controlcard,"MAXMIN",maxmin,100) call readi(controlcard,"MAXFUN",maxfun,100) call reada(controlcard,"TOLF",tolf,1.0d-6) call reada(controlcard,"RTOLF",rtolf,1.0d-6) out_minim=0 if (index(controlcard,"OUT_MINIM").gt.0) out_minim=iout print_ini=0 if (index(controlcard,"PRINT_INI").gt.0) print_ini=1 print_fin=0 if (index(controlcard,"PRINT_FIN").gt.0) print_fin=1 print_stat=0 if (index(controlcard,"PRINT_STAT").gt.0) print_stat=1 call reada(controlcard,"RSTIME",rstime,9.0d2) call reads(controlcard,"MINIMIZER",minimizer,"SUMSL") call readi(controlcard,"OPT_MODE",opt_mode,0) mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0 if (read_stat.eq.0 .and. mod_other_params) then write (iout,*) "Error: only optimization of energy-term", & " weights can be performed with READ_STAT=",read_stat call flush(iout) return1 endif if (index(controlcard,"RESTART").gt.0) then restart_flag=.true. else restart_flag=.false. endif c print *,me," checkpoint 6" return end c------------------------------------------------------------------------------------------------- subroutine read_pmf_data(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,kolor,klucz #endif include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.VMCPAR" include "COMMON.TORSION" include "COMMON.INTERACT" include "COMMON.ENERGIES" include "COMMON.MINPAR" include "COMMON.IOUNITS" include "COMMON.TIME1" include "COMMON.PROTFILES" include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.THERMAL" include "COMMON.FFIELD" include "COMMON.OPTIM" include "COMMON.CONTROL" include "COMMON.SCCOR" include "COMMON.SPLITELE" character*800 controlcard integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le integer ind,itype1,itype2,itypf,itypsc,itypp,itypt1,itypt2 integer ilen,lenpot,lenpre external ilen character*4 liczba,liczba1 #ifndef ISNAN external proc_proc #endif character*16 ucase character*16 key external ucase double precision xchuj,weitemp,weitemp_low,weitemp_up character*80 item(7) integer nitem integer rescode #ifdef NEWCORR c A 2/11/18 Read PMF parameters call card_concat(controlcard,.true.) torsion_pmf=index(controlcard,"TORSION_PMF").gt.0 turn_pmf=index(controlcard,"TURN_PMF").gt.0 eello_pmf=index(controlcard,"EELLO_PMF").gt.0 write (iout,*) "TORSION_PMF", TORSION_PMF," TURN_PMF ",TURN_PMF, & " EELLO_PMF",EELLO_PMF call reada(controlcard,"WELLO_PMF",wello_pmf,1.0d0) call reada(controlcard,"WELLO_PMF_LOW",wello_pmf_low,0.0d0) call reada(controlcard,"WELLO_PMF_UP",wello_pmf_up,5.0d0) call reada(controlcard,"WTURN_PMF",wturn_pmf,1.0d0) call reada(controlcard,"WTURN_PMF_LOW",wturn_pmf_low,0.0d0) call reada(controlcard,"WTURN_PMF_UP",wturn_pmf_up,5.0d0) call reada(controlcard,"WPMF",wpmf,1.0d0) write (iout,*) "nloctyp",nloctyp call multreada(controlcard,"E0",e0(0,-nloctyp+1), & (ntyp+1)*(2*nloctyp-1),0.0d0) call multreada(controlcard,"E0_LOW",e0_low(0,-nloctyp+1), & (ntyp+1)*(2*nloctyp-1), & -1.0d2) call multreada(controlcard,"E0_UP",e0_up(0,-nloctyp+1), & (ntyp+1)*(2*nloctyp-1), & 1.0d2) write (iout,*) "WTURN_PMF",WTURN_PMF,WTURN_PMF_LOW,WTURN_PMF_UP write (iout,*) "WELLO_PMF",WELLO_PMF,WELLO_PMF_LOW,WELLO_PMF_UP write (iout,*) "E0" do i=0,2 do j=-2,2 write (iout,"(2i5,3f15.5)") i,j,e0(i,j),e0_low(i,j),e0_up(i,j) enddo enddo #endif return end c----------------------------------------------------------------------- subroutine read_optim_parm(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,kolor,klucz #endif include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.VMCPAR" include "COMMON.TORSION" include "COMMON.LOCAL" include "COMMON.INTERACT" include "COMMON.ENERGIES" include "COMMON.MINPAR" include "COMMON.IOUNITS" include "COMMON.TIME1" include "COMMON.PROTFILES" include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.THERMAL" include "COMMON.FFIELD" include "COMMON.OPTIM" include "COMMON.CONTROL" include "COMMON.SCCOR" character*800 controlcard character*4 reskind integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le integer ind,ind1,ind2,itype1,itype2,itypf,itypsc,itypp, & itypt1,itypt2,masktemp,iblock,iaux,itypa integer ilen,lenpot,lenpre external ilen double precision aux,vb_low,vb_up,vb_rel character*4 liczba,liczba1 #ifndef ISNAN external proc_proc #endif character*16 ucase character*16 key external ucase double precision xchuj,weitemp,weitemp_low,weitemp_up character*80 item(7) character*3 typf,typa integer nitem integer rescode integer ind_shield /25/ lenpot=ilen(pot) lenpre=ilen(prefix) write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params do iblock=1,2 do itypt1=-ntyp,ntyp do itypt2=-ntyp,ntyp mask_tor(0,itypt1,itypt2,iblock)=0 weitor(0,itypt1,itypt2,iblock)=1.0d0 weitor_low(0,itypt1,itypt2,iblock)=1.0d0 weitor_up(0,itypt1,itypt2,iblock)=1.0d0 enddo enddo do itypt1=-ntyp,ntyp do itypt2=-ntyp,ntyp do l=0,3 mask_tor(l,itypt1,itypt2,iblock)=0 weitor(l,itypt1,itypt2,iblock)=1.0 weitor_low(l,itypt1,itypt2,iblock)=1.0 weitor_up(l,itypt1,itypt2,iblock)=1.0 enddo enddo enddo enddo #ifdef DEBUG write (iout,*) "ntyp",ntyp do itypt1=-ntyp,ntyp do itypt2=-ntyp,ntyp write (iout,*) "itypt1",itypt1," itypt2",itypt2, & "weitor",weitor(0,itypt1,itypt2,1),eitor(0,itypt1,itypt2,2) enddo enddo #endif if (mod_other_params) then c write (*,*) c &"Internal parameters of UNRES energy components will be optimized" call card_concat(controlcard,.true.) write (iout,*) "mod_side ",controlcard call flush(iout) mod_side = (index(controlcard,"MOD_SIDE").gt.0) if (mod_side) then nsingle_sc=0 npair_sc=0 call card_concat(controlcard,.true.) do while ( index(controlcard,"END").eq.0 ) call split_string(controlcard,item(1),4,nitem) if (nitem.eq.1 .or. item(2)(:1).eq."*") then nsingle_sc=nsingle_sc+1 ityp_ssc(nsingle_sc)=rescode(1,item(1),0,.false.) if (nitem.lt.3) then epss_low(nsingle_sc)=0.02d0 else read (item(3),*) epss_low(nsingle_sc) endif if (nitem.lt.4) then epss_up(nsingle_sc)=0.0d0 else read (item(4),*) epss_up(nsingle_sc) endif else npair_sc=npair_sc+1 ityp_psc(1,npair_sc)=rescode(1,item(1),0,.false.) ityp_psc(2,npair_sc)=rescode(1,item(2),0,.false.) if (nitem.lt.3) then epsp_low(npair_sc)=0.02d0 else read (item(3),*) epsp_low(npair_sc) endif if (nitem.lt.4) then epsp_up(npair_sc)=0.0d0 else read (item(4),*) epsp_up(npair_sc) endif endif call card_concat(controlcard,.true.) enddo if (nsingle_sc+npair_sc.eq.0) mod_side=.false. call card_concat(controlcard,.true.) endif mod_side_other= & index(controlcard,"MOD_SIDE_OTHER").gt.0 write (iout,*) "mod_side_other ",controlcard,mod_side_other if (mod_side_other) then mod_side_other=.false. do i=1,ntyp do j=1,5 mask_sigma(j,i)=0 enddo enddo call card_concat(controlcard,.true.) c write (iout,*) "mod_side_oter ",controlcard do while ( index(controlcard,"END").eq.0 ) call reads(controlcard,"RESKIND",reskind," ") itypsc=rescode(1,reskind,0,.false.) if (itypsc.lt.1 .or. itypsc.gt.20) then write (iout,*) "Error in SC optimization data;", & " SC type must not be dummy, type is" ,restyp," ",itypsc write (*,*) "Error in SC optimization data;", & " SC type must not be dummy, type is" ,restyp," ",itypsc return1 endif call readi(controlcard,"MASK_SIGMA0",mask_sigma(1,itypsc),0) call readi(controlcard,"MASK_SIGMAII",mask_sigma(2,itypsc),0) call readi(controlcard,"MASK_CHIP",mask_sigma(3,itypsc),0) call readi(controlcard,"MASK_ALP",mask_sigma(4,itypsc),0) call reada(controlcard,"SIGMA0_LOW",sigma_low(1,itypsc),0.d0) call reada(controlcard,"SIGMA0_UP",sigma_up(1,itypsc),0.0d0) call reada(controlcard,"SIGMAII_LOW",sigma_low(2,itypsc), & 0.0d0) call reada(controlcard,"SIGMAII_UP",sigma_up(2,itypsc),0.0d0) call reada(controlcard,"CHIP_LOW",sigma_low(3,itypsc),0.d0) call reada(controlcard,"CHIP_UP",sigma_up(3,itypsc),0.0d0) call reada(controlcard,"ALP_LOW",sigma_low(4,itypsc),0.d0) call reada(controlcard,"ALP_UP",sigma_up(4,itypsc),0.0d0) do k=1,4 if (sigma_low(k,itypsc).eq.0.0d0 .and. & sigma_up(k,itypsc).eq.0.0d0) mask_sigma(k,itypsc)=0 enddo do k=1,4 mod_side_other=mod_side_other.or.mask_sigma(k,itypsc).gt.0 enddo write (iout,'(a4,i3,4(i2,2f8.3))') reskind,itypsc, & (mask_sigma(k,itypsc), & sigma_low(k,itypsc),sigma_up(k,itypsc),k=1,4) call card_concat(controlcard,.true.) c write (iout,*) "mod_side_oter ",controlcard enddo write (iout,*) "Optimization flags of one-body SC parameters" do i=1,ntyp write (iout,'(a4,i3,4(i2,2f8.3))') restyp(i),i, & (mask_sigma(k,i),sigma_low(k,i),sigma_up(k,i),k=1,4) enddo call card_concat(controlcard,.true.) endif c write (iout,*) "mod_side_other ",mod_side_other c write (iout,*) "mod_tor ",controlcard c call flush(iout) mod_tor = index(controlcard,"MOD_TOR").gt.0 if (mod_tor) then do iblock=1,2 do i=-ntortyp,ntortyp do j=-ntortyp,ntortyp mask_tor(0,i,j,iblock)=0 weitor(0,i,j,iblock)=1.0d0 weitor_low(0,i,j,iblock)=0.0d0 weitor_up(0,i,j,iblock)=2.0d0 enddo enddo enddo call card_concat(controlcard,.true.) write (iout,*) controlcard do while ( index(controlcard,"END").eq.0 ) call split_string(controlcard,item(1),7,nitem) if (nitem.lt.2) then write (*,*) "Error in torsional optimization data;", & " must specify both residues and type" return1 endif weitemp=1.0d0 weitemp_low=0.0d0 weitemp_up=2.0d0 masktemp=1 iblock=1 write (iout,*) "item3 ",item(3)," item4 ",item(4), & " item5",item(5) call flush(iout) if (nitem.gt.2) read(item(3),*) masktemp if (nitem.gt.3) read(item(4),*) weitemp if (nitem.gt.4) read(item(5),*) weitemp_low if (nitem.gt.5) read(item(6),*) weitemp_up if (nitem.gt.6) read(item(7),*) iblock write (iout,*) controlcard write (iout,*) item(1)," ",item(2),weitemp, & weitemp_low,weitemp_up if (index(item(1),"*").gt.0) then ist1=1 iet1=ntyp else ist1=itortyp(rescode(1,item(1),0,.false.)) iet1=ist1 endif if (index(item(2),"*").gt.0) then ist2=1 iet2=ntyp else ist2=itortyp(rescode(1,item(2),0,.false.)) iet2=ist2 endif c write (iout,*) "ist1",ist1," iet1",iet1," ist2",ist2, c & " iet2",iet2 c call flush(iout) do itypt1=ist1,iet1 do itypt2=ist2,iet2 c write (iout,*) "itypt1",itypt1," itypt2",itypt2, c & weitemp,weitemp_low,weitemp_up mask_tor(0,itypt1,itypt2,iblock)=masktemp weitor(0,itypt1,itypt2,iblock)=weitemp weitor_low(0,itypt1,itypt2,iblock)=weitemp_low weitor_up(0,itypt1,itypt2,iblock)=weitemp_up c write (iout,*) "itypt1",itypt1," itypt2",itypt2, c & mask_tor(0,itypt1,itypt2,iblock), c & weitor(0,itypt1,itypt2,iblock), c & weitor_low(0,itypt1,itypt2,iblock), c & weitor_up(0,itypt1,itypt2,iblock) enddo enddo call card_concat(controlcard,.true.) write (iout,*) controlcard enddo #ifdef NEWCORR if (tor_mode.gt.1) then write (iout,*) "TORMODE is",tor_mode, & " torsional constants are computed from energy map expansion." endif #endif #ifdef DEBUG write (iout,*) "Initialized torsional parameters:" do iblock=1,2 do itypt1=-ntortyp,ntortyp do itypt2=-ntortyp,ntortyp write (iout,'(3i5,3f10.5)') itypt1,itypt2, & mask_tor(0,itypt1,itypt2,iblock), & weitor(0,itypt1,itypt2,iblock), & weitor_low(0,itypt1,itypt2,iblock), & weitor_up(0,itypt1,itypt2,iblock) enddo enddo enddo #endif if (tor_mode.eq.1) then c Exchange indices because the residue order in new torsionals is reversed do iblock=1,2 do itypt1=-ntortyp,ntortyp do itypt2=itypt1+1,ntortyp iaux=mask_tor(0,itypt1,itypt2,iblock) mask_tor(0,itypt1,itypt2,iblock)= & mask_tor(0,itypt2,itypt1,iblock) mask_tor(0,itypt2,itypt1,iblock)=iaux aux=weitor(0,itypt1,itypt2,iblock) weitor(0,itypt1,itypt2,iblock)= & weitor(0,itypt2,itypt1,iblock) weitor(0,itypt2,itypt1,iblock)=aux aux=weitor_low(0,itypt1,itypt2,iblock) weitor_low(0,itypt1,itypt2,iblock)= & weitor_low(0,itypt2,itypt1,iblock) weitor_low(0,itypt2,itypt1,iblock)=aux aux=weitor_up(0,itypt1,itypt2,iblock) weitor_up(0,itypt1,itypt2,iblock)= & weitor_up(0,itypt2,itypt1,iblock) weitor_up(0,itypt2,itypt1,iblock)=aux enddo enddo enddo endif call card_concat(controlcard,.true.) endif write (iout,*) "mod_sccor ",controlcard call flush(iout) mod_sccor = index(controlcard,"MOD_SCCOR").gt.0 if (mod_sccor) then call card_concat(controlcard,.true.) do iblock=1,2 do l=1,3 do i=-nsccortyp,nsccortyp do j=-nsccortyp,nsccortyp mask_tor(l,i,j,iblock)=0 weitor(l,i,j,iblock)=1.0d0 weitor_low(l,i,j,iblock)=0.0d0 weitor_up(l,i,j,iblock)=2.0d0 enddo enddo enddo enddo do while ( index(controlcard,"END").eq.0 ) call split_string(controlcard,item(1),7,nitem) if (nitem.lt.3) then write (*,*) "Error in torsional optimization data;", & " must specify both residues and type" return1 endif weitemp=1.0d0 weitemp_low=0.0d0 weitemp_up=0.0d0 if (nitem.gt.3) read(item(4),*) masktemp if (nitem.gt.4) read(item(5),*) weitemp if (nitem.gt.5) read(item(6),*) weitemp_low if (nitem.gt.6) read(item(7),*) weitemp_up if (index(item(1),"*").gt.0) then ist1=1 iet1=ntyp else ist1=isccortyp(rescode(1,item(1),0,.false.)) iet1=ist1 endif if (index(item(2),"*").gt.0) then ist2=1 iet2=ntyp else ist2=isccortyp(rescode(1,item(2),0,.false.)) iet2=ist2 endif if (index(item(3),"*").gt.0) then ls=1 le=3 else read(item(3),*) ls le=ls endif do itypt1=ist1,iet1 do itypt2=ist2,iet2 do l=ls,le mask_tor(l,itypt1,itypt2,1)=masktemp weitor(l,itypt1,itypt2,1)=weitemp weitor_low(l,itypt1,itypt2,1)=weitemp_low weitor_up(l,itypt1,itypt2,1)=weitemp_up enddo enddo enddo call card_concat(controlcard,.true.) enddo call card_concat(controlcard,.true.) endif #ifdef DEBUG write (iout,*) "Optimizable sidechain-torsional parameters:" do itypt1=-nsccortyp,nsccortyp do itypt2=-nsccortyp,nsccortyp do l=1,3 if (mask_tor(l,itypt1,itypt2,1).gt.0) & write (iout,'(4i5,3f10.5)') itypt1,itypt2,l, & mask_tor(l,itypt1,itypt2,1),weitor(l,itypt1,itypt2,1), & weitor_low(l,itypt1,itypt2,1),weitor_up(l,itypt1,itypt2,1) enddo enddo enddo #endif mod_ang=tor_mode.gt.0 .and. index(controlcard,"MOD_ANGLE").gt.0 write (iout,*) "mod_angle ",controlcard write (iout,*) "mod_ang",mod_ang if (mod_ang) then do i=0,nthetyp-1 mask_ang(i)=0 enddo call card_concat(controlcard,.true.) do while (index(controlcard,'END').eq.0) write (iout,'(a)') "angle: ",controlcard call reads(controlcard,"TYPE",typa," ") itypa=rescode(1,typa,0,.false.) c write (iout,*) "TYPA ",typa," itypa",itypa if (itypa.eq.0 .or. itypa.gt.ntyp) then write (*,*) "Error specifying angle parms" stop endif itypa=itortyp(itypa) mask_ang(itypa)=1 write (iout,*) "bend type",itypa call reada(controlcard,"VB_LOW",vb_low,-1.0d5) call reada(controlcard,"VB_UP",vb_up,1.0d5) call reada(controlcard,"VB_REL",vb_rel,0.0d0) write (iout,*) "VB_LOW",vb_low," VB_UP",vb_up," VB_REL",vb_rel do i=1,nbend_kcc_TB(itypa) if (vb_rel.gt.0) then write (iout,*) "v1bend_chyb=",v1bend_chyb(i,itypa) v1bend_low(i,itypa)=v1bend_chyb(i,itypa)*(1.0d0 & -dsign(vb_rel,v1bend_chyb(i,itypa))) v1bend_up(i,itypa)=v1bend_chyb(i,itypa)*(1.0d0 & +dsign(vb_rel,v1bend_chyb(i,itypa))) else v1bend_low(i,itypa)=vb_low v1bend_up(i,itypa)=vb_up endif enddo call card_concat(controlcard,.true.) enddo call card_concat(controlcard,.true.) write (iout,*) "Boundaries of angle potential coefficients" write (iout,'(3a)') "Index"," Low bound"," Up bound" do i=0,nthetyp-1 if (mask_ang(i).eq.0) cycle write (iout,'(2a)') 'Type ',restyp(iloctyp(i)) do j=1,nbend_kcc_TB(i) write (iout,'(i5,2f15.1)') j,v1bend_low(j,i),v1bend_up(j,i) enddo enddo endif write (iout,*) "mod_fourier ",controlcard call flush(iout) mod_fourier(nloctyp)=index(controlcard,"MOD_FOURIER") #ifdef NEWCORR if (mod_fourier(nloctyp).gt.0) then mod_fourier(nloctyp)=0 do i=0,nloctyp-1 mod_fourier(i)=0 do ii=1,3 do j=1,2 mask_bnew1(ii,j,i)=0 mask_bnew2(ii,j,i)=0 mask_ccnew(ii,j,i)=0 mask_ddnew(ii,j,i)=0 enddo enddo do ii=1,2 do j=1,2 do k=1,2 mask_eenew(ii,j,k,i)=0 enddo enddo enddo do ii=1,3 mask_e0new(ii,i)=0 enddo enddo call card_concat(controlcard,.true.) do while ( index(controlcard,"END").eq.0 ) c write(iout,*) controlcard call reads(controlcard,"TYPF",typf," ") itypf=rescode(1,typf,0,.false.) c write (iout,*) "TYPF ",typf," itypf",itypf if (itypf.eq.0 .or. itypf.gt.ntyp) then write (*,*) "Error specifying fourier parms" stop endif itypf=itype2loc(itypf) write (iout,*) "local type",itypf if (index(controlcard,"B1_LOW").gt.0) then call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying B1, components undefined",ind,ii stop endif mask_bnew1(ii,ind,itypf)=1 call reada(controlcard,"B1_LOW",bnew1_low(ii,ind,itypf), & 0.1d0) call reada(controlcard,"B1_UP",bnew1_up(ii,ind,itypf), & 0.0d0) mod_fourier(itypf)=mod_fourier(itypf) & +mask_bnew1(ii,ind,itypf) else if (index(controlcard,"B2_LOW").gt.0) then mask_bnew2(ii,ind,itypf)=1 call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying B2, components undefined",ind,ii stop endif call reada(controlcard,"B2_LOW",bnew2_low(ii,ind,itypf), & 0.1d0) call reada(controlcard,"B2_UP",bnew2_up(ii,ind,itypf), & 0.0d0) mod_fourier(itypf)=mod_fourier(itypf) & +mask_bnew2(ii,ind,itypf) else if (index(controlcard,"C_LOW").gt.0) then call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying C, components undefined",ind,ii stop endif mask_ccnew(ii,ind,itypf)=1 call reada(controlcard,"C_LOW",ccnew_low(ii,ind,itypf),.1d0) call reada(controlcard,"C_UP",ccnew_up(ii,ind,itypf),0.0d0) mod_fourier(itypf)=mod_fourier(itypf) & +mask_ccnew(ii,ind,itypf) else if (index(controlcard,"D_LOW").gt.0) then call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying D, components undefined",ind,ii stop endif mask_ddnew(ii,ind,itypf)=1 call reada(controlcard,"D_LOW",ddnew_low(ii,ind,itypf), & 0.1d0) call reada(controlcard,"D_UP",ddnew_up(ii,ind,itypf), & 0.0d0) mod_fourier(itypf)=mod_fourier(itypf) & +mask_ddnew(ii,ind,itypf) else if (index(controlcard,"E0_LOW").gt.0) then call readi(controlcard,"COEFF",ii,0) if (ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying E0, components undefined",ii stop endif mask_e0new(ii,itypf)=1 call reada(controlcard,"E0_LOW",e0new_low(ii,itypf), & 0.1d0) call reada(controlcard,"E0_UP",e0new_up(ii,itypf), & 0.0d0) mod_fourier(itypf)=mod_fourier(itypf) & +mask_e0new(ii,itypf) else if (index(controlcard,"E_LOW").gt.0) then call readi(controlcard,"IND1",ind1,0) call readi(controlcard,"IND2",ind2,0) call readi(controlcard,"COEFF",ii,0) if (ind1.eq.0 .or. ind2.eq.0 .or. ii.eq.0 .or. ii.gt.2) then write (iout,*) & "Error specifying E, components undefined",ind1,ind2,ii stop endif mask_eenew(ii,ind1,ind2,itypf)=1 call reada(controlcard,"E_LOW", & eenew_low(ii,ind1,ind2,itypf),0.1d0) call reada(controlcard,"E_UP", & eenew_up(ii,ind1,ind2,itypf),0.0d0) mod_fourier(itypf)=mod_fourier(itypf) & +mask_eenew(ii,ind1,ind2,itypf) endif call card_concat(controlcard,.true.) enddo call card_concat(controlcard,.true.) write (iout,*) "mod_fouriertor card ",controlcard mod_fouriertor(nloctyp)=index(controlcard,"MOD_FOURIERTOR") write (iout,*) "mod_fouriertor value",mod_fouriertor(nloctyp) write (2,*) "SPLIT_FOURIERTOR",SPLIT_FOURIERTOR, & " tor_mode",tor_mode IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2 & .and. mod_fouriertor(nloctyp).gt.0) THEN do i=0,nloctyp-1 mod_fouriertor(i)=0 do ii=1,3 do j=1,2 mask_bnew1tor(ii,j,i)=0 mask_bnew2tor(ii,j,i)=0 mask_ccnewtor(ii,j,i)=0 mask_ddnewtor(ii,j,i)=0 enddo enddo do ii=1,2 do j=1,2 do k=1,2 mask_eenewtor(ii,j,k,i)=0 enddo enddo enddo do ii=1,3 mask_e0newtor(ii,i)=0 enddo enddo call card_concat(controlcard,.true.) do while ( index(controlcard,"END").eq.0 ) c write(iout,*) controlcard call reads(controlcard,"TYPF",typf," ") itypf=rescode(1,typf,0,.false.) c write (iout,*) "TYPF ",typf," itypf",itypf if (itypf.eq.0 .or. itypf.gt.ntyp) then write (*,*) "Error specifying fourier parms" stop endif itypf=itype2loc(itypf) write (iout,*) "local type",itypf if (index(controlcard,"B1_LOW").gt.0) then call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying B1, components undefined",ind,ii stop endif mask_bnew1tor(ii,ind,itypf)=1 call reada(controlcard,"B1_LOW",bnew1tor_low(ii,ind,itypf), & 0.1d0) call reada(controlcard,"B1_UP",bnew1tor_up(ii,ind,itypf), & 0.0d0) mod_fouriertor(itypf)=mod_fouriertor(itypf) & +mask_bnew1tor(ii,ind,itypf) else if (index(controlcard,"B2_LOW").gt.0) then mask_bnew2tor(ii,ind,itypf)=1 call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying B2, components undefined",ind,ii stop endif call reada(controlcard,"B2_LOW",bnew2tor_low(ii,ind,itypf), & 0.1d0) call reada(controlcard,"B2_UP",bnew2tor_up(ii,ind,itypf), & 0.0d0) mod_fouriertor(itypf)=mod_fouriertor(itypf) & +mask_bnew2tor(ii,ind,itypf) else if (index(controlcard,"C_LOW").gt.0) then call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying C, components undefined",ind,ii stop endif mask_ccnewtor(ii,ind,itypf)=1 call reada(controlcard,"C_LOW",ccnewtor_low(ii,ind,itypf), & .1d0) call reada(controlcard,"C_UP",ccnewtor_up(ii,ind,itypf), & 0.0d0) mod_fouriertor(itypf)=mod_fouriertor(itypf) & +mask_ccnewtor(ii,ind,itypf) else if (index(controlcard,"D_LOW").gt.0) then call readi(controlcard,"IND",ind,0) call readi(controlcard,"COEFF",ii,0) if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying D, components undefined",ind,ii stop endif mask_ddnewtor(ii,ind,itypf)=1 call reada(controlcard,"D_LOW",ddnewtor_low(ii,ind,itypf), & 0.1d0) call reada(controlcard,"D_UP",ddnewtor_up(ii,ind,itypf), & 0.0d0) mod_fouriertor(itypf)=mod_fouriertor(itypf) & +mask_ddnewtor(ii,ind,itypf) else if (index(controlcard,"E0_LOW").gt.0) then call readi(controlcard,"COEFF",ii,0) if (ii.eq.0 .or. ii.gt.3) then write (iout,*) & "Error specifying E0, components undefined",ii stop endif mask_e0newtor(ii,itypf)=1 call reada(controlcard,"E0_LOW",e0newtor_low(ii,itypf), & 0.1d0) call reada(controlcard,"E0_UP",e0newtor_up(ii,itypf), & 0.0d0) mod_fouriertor(itypf)=mod_fouriertor(itypf) & +mask_e0newtor(ii,itypf) else if (index(controlcard,"E_LOW").gt.0) then call readi(controlcard,"IND1",ind1,0) call readi(controlcard,"IND2",ind2,0) call readi(controlcard,"COEFF",ii,0) if (ind1.eq.0 .or. ind2.eq.0 .or. ii.eq.0 .or. ii.gt.2) then write (iout,*) & "Error specifying E, components undefined",ind1,ind2,ii stop endif mask_eenewtor(ii,ind1,ind2,itypf)=1 call reada(controlcard,"E_LOW", & eenewtor_low(ii,ind1,ind2,itypf),0.1d0) call reada(controlcard,"E_UP", & eenewtor_up(ii,ind1,ind2,itypf),0.0d0) mod_fouriertor(itypf)=mod_fouriertor(itypf) & +mask_eenewtor(ii,ind1,ind2,itypf) endif call card_concat(controlcard,.true.) enddo call card_concat(controlcard,.true.) write (2,*) "End read FOURIERTOR ",controlcard ENDIF ! SPLIT_FOURIERTOR endif do itypf=0,nloctyp-1 write (iout,*) "itypf",itypf," mod_fourier", & mod_fourier(itypf) mod_fourier(nloctyp)=mod_fourier(nloctyp) & +mod_fourier(itypf) enddo write (iout,*) "mod_fourier all",mod_fourier(nloctyp) do itypf=0,nloctyp-1 write (iout,*) "itypf",itypf," mod_fouriertor", & mod_fouriertor(itypf) mod_fouriertor(nloctyp)=mod_fouriertor(nloctyp) & +mod_fouriertor(itypf) enddo write (iout,*) "mod_fouriertor all",mod_fouriertor(nloctyp) #else if (mod_fourier(nloctyp).gt.0) then call card_concat(controlcard,.true.) do while ( index(controlcard,"END").eq.0 ) call readi(controlcard,"TYPF",typf," ") itypf=rescode(1,typf,0,.false.) if (itypf.eq.0 .or. itypf.gt.ntyp) then write (*,*) "Error specifying fourier parms" stop endif itypf=itype2loc(itypf) call readi(controlcard,"IND",ind,0) call reada(controlcard,"B_LOW",b_low(ind,itypf),0.1d0) call reada(controlcard,"B_UP",b_up(ind,itypf),0.0d0) mask_fourier(ind,itypf)=1 mod_fourier(itypf)=mod_fourier(itypf) & +mask_fourier(ind,itypf) mod_fourier(nloctyp)=mod_fourier(nloctyp) & +mask_fourier(ind,itypf) call card_concat(controlcard,.true.) enddo call card_concat(controlcard,.true.) endif do itypf=0,nloctyp-1 write (iout,*) "itypf",itypf," mod_fourier",mod_fourier(itypf) mod_fourier(nloctyp)=mod_fourier(nloctyp)+mod_fourier(itypf) enddo write (iout,*) "mod_fourier all",mod_fourier(nloctyp) #endif do i=1,2 do j=1,2 do k=1,2 mask_elec(i,j,k)=0 enddo enddo enddo mod_elec=index(controlcard,"MOD_ELEC").gt.0 if (mod_elec) then mod_elec=.false. call card_concat(controlcard,.true.) do while ( index(controlcard,"END").eq.0 ) c write (iout,*) "controlcard ",controlcard call readi(controlcard,"TYPE1",itype1,0) call readi(controlcard,"TYPE2",itype2,0) c write (iout,*) "itype1",itype1," itype2",itype2 if (itype1.eq.0 .or. itype1.gt.2 .or. itype2.eq.0 & .or. itype2.gt.2) then write (iout,*) "Error specifying elec parms" stop endif if (index(controlcard,"EPP").gt.0) then mod_elec=.true. mask_elec(itype1,itype2,1)=1 mask_elec(itype2,itype1,1)=1 call reada(controlcard,"EPP_LOW",epp_low(itype1,itype2), & 0.1d0) epp_low(itype2,itype1)=epp_low(itype1,itype2) call reada(controlcard,"EPP_UP",epp_up(itype1,itype2), & 0.0d0) epp_up(itype2,itype1)=epp_up(itype1,itype2) endif if (index(controlcard,"RPP").gt.0) then mod_elec=.true. mask_elec(itype1,itype2,2)=1 mask_elec(itype2,itype1,2)=1 call reada(controlcard,"RPP_LOW",rpp_low(itype1,itype2), & 0.1d0) rpp_low(itype2,itype1)=rpp_low(itype1,itype2) call reada(controlcard,"RPP_UP",rpp_up(itype1,itype2), & 0.0d0) rpp_up(itype2,itype1)=rpp_up(itype1,itype2) endif if (index(controlcard,"ELPP6").gt.0) then mod_elec=.true. mask_elec(itype1,itype2,3)=1 mask_elec(itype2,itype1,3)=1 call reada(controlcard,"ELPP6_LOW", & elpp6_low(itype1,itype2),0.1d0) elpp6_low(itype2,itype1)=elpp6_low(itype1,itype2) call reada(controlcard,"ELPP6_UP", & elpp6_up(itype1,itype2),0.0d0) elpp6_up(itype2,itype1)=elpp6_up(itype1,itype2) endif if (index(controlcard,"ELPP3").gt.0) then mod_elec=.true. mask_elec(itype1,itype2,4)=1 mask_elec(itype2,itype1,4)=1 call reada(controlcard,"ELPP3_LOW", & elpp3_low(itype1,itype2),0.1d0) elpp3_low(itype2,itype1)=elpp3_low(itype1,itype2) call reada(controlcard,"ELPP3_UP", & elpp3_up(itype1,itype2),0.0d0) elpp3_up(itype2,itype1)=elpp3_up(itype1,itype2) endif call card_concat(controlcard,.true.) enddo call card_concat(controlcard,.true.) endif do i=1,20 do j=1,2 do k=1,2 mask_scp(i,j,k)=0 enddo enddo enddo mod_scp=index(controlcard,"MOD_SCP").gt.0 if (mod_scp) then mod_scp=.false. call card_concat(controlcard,.true.) do while (index(controlcard,"END").eq.0) if (index(controlcard,"EPSCP").gt.0) then mod_scp=.true. call readi(controlcard,"ITYPSC",itypsc,0) call readi(controlcard,"ITYPP",itypp,0) if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then write (iout,*) "Error specifying scp parms" stop endif mask_scp(itypsc,itypp,1)=1 call reada(controlcard,"EPSCP_LOW", & epscp_low(itypsc,itypp),0.1d0) call reada(controlcard,"EPSCP_UP", & epscp_up(itypsc,itypp),0.0d0) endif if (index(controlcard,"RSCP").gt.0) then mod_scp=.true. call readi(controlcard,"ITYPSC",itypsc,0) call readi(controlcard,"ITYPP",itypp,0) mask_scp(itypsc,itypp,2)=1 call readi(controlcard,"ITYPSC",itypsc,0) call readi(controlcard,"ITYPP",itypp,0) if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then write (iout,*) "Error specifying scp parms" stop endif call reada(controlcard,"RSCP_LOW", & rscp_low(itypsc,itypp),0.1d0) c write(iout,*)itypsc,itypp,rscp_low(itypsc,itypp) call reada(controlcard,"RSCP_UP", & rscp_up(itypsc,itypp),0.0d0) c write(iout,*)itypsc,itypp,rscp_up(itypsc,itypp) endif call card_concat(controlcard,.true.) enddo endif write (iout,*) "END ",controlcard call flush(iout) mod_other_params= & mod_fourier(nloctyp).gt.0 .or. mod_elec .or. mod_scp & .or. mod_sccor .or. mod_ang .or. imask(ind_shield).eq.1 if (read_stat.lt.2. .and. mod_other_params) then write (iout,*)"Error - only weights and sidechain parameters", & " can be optimized with READ_STAT=",read_stat call flush(iout) return1 endif init_ene = init_ene .or. read_stat.eq.2 endif write (iout,*) "End read: MOD_OTHER_PARAMS ",mod_other_params c write (*,*) "MOD_SIDE ",mod_side," MOD_FOURIER", c & mod_fourier(nloctyp), c & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp init_ene=init_ene .or. mod_other_params c write (iout,*) "init_ene",init_ene c write (iout,*) call flush(iout) return end c----------------------------------------------------------------------------- subroutine print_general_data(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,kolor,klucz #endif include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.VMCPAR" include "COMMON.TORSION" include "COMMON.LOCAL" include "COMMON.INTERACT" include "COMMON.ENERGIES" include "COMMON.MINPAR" include "COMMON.IOUNITS" include "COMMON.TIME1" include "COMMON.PROTFILES" include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.THERMAL" include "COMMON.FFIELD" include "COMMON.OPTIM" include "COMMON.CONTROL" include "COMMON.SCCOR" character*800 controlcard integer i,j,k,l,ii,n_ene_found,itypt,jtypt integer ind,itype1,itype2,itypf,itypsc,itypp integer ilen,lenpot,lenpre external ilen character*4 liczba,liczba1 if (mode.eq.1) then write (iout,*) "Single-point target function evaluation" else if (mode.eq.2) then write (iout,*) "Grid search of the parameter space" else if (mode.eq.3) then write (iout,*) "Minimization of the target function" else write (iout,*) "Wrong MODE",mode,", should be 1, 2, or 3" call flush(iout) return1 endif write (iout,*) "RESCALE_MODE",rescale_mode c mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0 c if (read_stat.eq.0 .and. mod_other_params) then c write (iout,*) "Error: only optimization of energy-term", c & " weights can be performed with READ_STAT=",read_stat c call flush(iout) c return1 c endif if (mode.eq.3) then write (iout,*) "Parameters of the SUMSL procedure:" write (iout,'(a,t15,i5)') "MAXMIN",maxmin write (iout,'(a,t15,i5)') "MAXFUN",maxfun write (iout,'(a,t15,e15.7)') "TOLF",tolf write (iout,'(a,t15,e15.7)') "RTOLF",rtolf endif if (mod_other_params) then write (iout,*) &"Internal parameters of UNRES energy components will be optimized" c call card_concat(controlcard,.true.) c mod_side = (index(controlcard,"MOD_SIDE").gt.0) if (mod_side) then write (iout,*) "SC epsilons to be optimized:" write (iout,*) "Single: eps(i,j)=eps(i,j)+(x(i)+x(j))/2" do i=1,nsingle_sc write (iout,*) restyp(ityp_ssc(i)),epss_low(i),epss_up(i) enddo write (iout,*) "Pairs: eps(i,j)=eps(i,j)+x(i,j):" do i=1,npair_sc write (iout,*) restyp(ityp_psc(1,i)), & restyp(ityp_psc(2,i)),epsp_low(i),epsp_up(i) enddo endif if (mod_sccor) then write (iout,*)"Torsional weights/coefficients to be optimized" write(iout,'(a)') "TYP1 TYP2 WEIGHT", & " LOWER BOUND UPPER BOUND" do itypt=-nsccortyp,nsccortyp do jtypt=-nsccortyp,nsccortyp do l=1,3 if (mask_tor(l,itypt,jtypt,1).gt.0) then write(iout,'(3a4,3f10.5)')l,restyp(itypt),restyp(jtypt), & weitor(l,itypt,jtypt,1),weitor_low(l,itypt,jtypt,1), & weitor_up(l,itypt,jtypt,1) endif enddo enddo enddo endif c 7/8/17 AL: Optimization of the bending parameters if (mod_ang) then write (iout,*) "Boundaries of angle potential coefficients" write (iout,'(3a)') "Index"," Low bound"," Up bound" do i=0,nthetyp if (mask_ang(i).eq.0) cycle write (iout,'(2a)') 'Type ',restyp(iloctyp(i)) do j=1,nbend_kcc_TB(i) write (iout,'(i5,2f15.3)') j,v1bend_low(j,i),v1bend_up(j,i) enddo enddo endif if (mod_fourier(nloctyp).gt.0) then write (iout,*) "Fourier coefficients to be optimized" do itypf=0,nloctyp-1 if (mod_fourier(itypf).gt.0) then write (iout,'(3a,i2)') "Type ",restyp(iloctyp(itypf)), & " number of coeffs",mod_fourier(itypf) write(iout,'(a)') "NAME COEFF LOWER BOUND UPPER BOUND" #ifdef NEWCORR do j=1,2 do k=1,3 if (mask_bnew1(k,j,itypf).gt.0) & write (iout,'(2hB1,2i2,f10.5)') k,j,bnew1(k,j,itypf), & bnew1_low(k,j,itypf),bnew1_up(k,j,itypf) enddo enddo do j=1,2 do k=1,3 if (mask_bnew2(k,j,itypf).gt.0) & write (iout,'(2hB2,2i2,3f11.5)') k,j,bnew2(k,j,itypf), & bnew2_low(k,j,itypf),bnew2_up(k,j,itypf) enddo enddo do j=1,2 do k=1,3 if (mask_ccnew(k,j,itypf).gt.0) & write (iout,'(2hCC,2i2,3f11.5)') k,j,ccnew(k,j,itypf), & ccnew_low(k,j,itypf),ccnew_up(k,j,itypf) enddo enddo do j=1,2 do k=1,3 if (mask_ddnew(k,j,itypf).gt.0) & write (iout,'(2hDD,2i2,3f11.5)') k,j,ddnew(k,j,itypf), & ddnew_low(k,j,itypf),ddnew_up(k,j,itypf) enddo enddo do j=1,2 if (mask_e0new(j,itypf).gt.0) & write (iout,'(2hE0,i2,3f11.5)') j,e0new(j,itypf), & e0new_low(j,itypf),e0new_up(j,itypf) enddo do j=1,2 do k=1,2 do l=1,2 if (mask_eenew(l,k,j,itypf).gt.0) & write (iout,'(2hEE,3i2,3f11.5)') & l,k,j,eenew(l,k,j,itypf),eenew_low(l,k,j,itypf), & eenew_up(l,k,j,itypf) enddo enddo enddo #else do i=2,13 if (mask_fourier(i,itypf).gt.0) then write (iout,'(1hB,i2,3f11.5)') & i,b(i,itypf),b_low(i,itypf),b_up(i,itypf) endif enddo #endif endif enddo endif if (mod_elec) then write (iout,*) write (iout,*) "Electrostatic parameters to be optimized" do itype1=1,2 do itype2=1,itype1 if (mask_elec(itype1,itype2,1).gt.0) & write (iout,'(2i3," EPP",f8.3," EPP_LOW",f8.3, & " EPP_UP",f8.3)') & itype1,itype2,epp(itype1,itype2),epp_low(itype1,itype2), & epp_up(itype1,itype2) if (mask_elec(itype1,itype2,2).gt.0) & write (iout,'(2i3," RPP",f8.3," RPP_LOW",f8.3, & " RPP_UP",f8.3)') & itype1,itype2,rpp(itype1,itype2),rpp_low(itype1,itype2), & rpp_up(itype1,itype2) if (mask_elec(itype1,itype2,3).gt.0) & write (iout,'(2i3," ELPP6",f8.3," ELPP6_LOW",f8.3, & " ELPP6_UP",f8.3)') & itype1,itype2,elpp6(itype1,itype2), & elpp6_low(itype1,itype2),elpp6_up(itype1,itype2) if (mask_elec(itype1,itype2,4).gt.0) & write (iout,'(2i3," ELPP3",f8.3," ELPP3_LOW",f8.3, & " ELPP3_UP",f8.3)') & itype1,itype2,elpp3(itype1,itype2), & elpp3_low(itype1,itype2),elpp3_up(itype1,itype2) enddo enddo endif if (mod_scp) then mod_scp=.false. write (iout,*) write (iout,*) "SCP parameters to be optimized:" if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+ & mask_scp(0,2,2) .gt. 0) then write (iout,*) & "SCP parameters are assumed to depend on peptide group type only" do j=1,2 if (mask_scp(0,j,1).gt.0) & write (iout,'(i3," EPSCP",f8.3," EPSCP_LOW",f8.3, & " EPSCP_UP",f8.3)') j,eps_scp(1,j),epscp_low(0,j), & epscp_up(0,j) if (mask_scp(0,j,2).gt.0) & write (iout,'(i3," RSCP",f8.3," RSCP_LOW",f8.3, & " RSCP_UP",f8.3)') j,rscp(1,j),rscp_low(0,j), & rscp_up(0,j) enddo else do i=1,20 do j=1,2 if (mask_scp(i,j,1).gt.0) & write (iout,'(2i3," EPSCP",f8.3," EPSCP_LOW",f8.3, & " EPSCP_UP",f8.3)') i,j,eps_scp(i,j),epscp_low(i,j), & epscp_up(i,j) if (mask_scp(i,j,2).gt.0) & write (iout,'(2i3," RSCP",f8.3," RSCP_LOW",f8.3, & " RSCP_UP",f8.3)') i,j,rscp(i,j),rscp_low(i,j), & rscp_up(i,j) enddo enddo endif endif endif write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params write (iout,*) "MOD_SIDE ",mod_side," MOD_FOURIER", & mod_fourier(nloctyp), & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp," mod_ang",mod_ang init_ene=init_ene .or. mod_other_params write (iout,*) "init_ene",init_ene write (iout,*) call flush(iout) return end c----------------------------------------------------------------------------- subroutine read_protein_data(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" #endif include "COMMON.CONTROL" include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.COMPAR" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.PROTNAME" include "COMMON.VMCPAR" include "COMMON.OPTIM" include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.ALLPROT" include "COMMON.THERMAL" include "COMMON.TIME1" include "COMMON.WEIGHTDER" character*64 nazwa,key character*16000 controlcard character*16000 all_protfiles integer i,j,k,kk,l,iprot,ii,if,ib,ibatch,nn,nn1,iww,maskcheck, & il,inat,ilevel,weightread,jfrag,jclass,nfragm,iparm integer nrec,nlines,iscor double precision energ,wangnorm(maxT),wq(maxT),wrms(maxT), & wrgy(maxT),wsign(maxT),wangnat(maxT),wqnat(maxT),wrmsnat(maxT), & wrgynat(maxT),wcorangnorm(maxT),wcorq(maxT), & wcorrms(maxT),wcorrgy(maxT),wcorsign(maxT),wcorangnat(maxT), & wcorqnat(maxT),wcorrmsnat(maxT),wcorrgynat(maxT), & angnormlow(maxT),qlow(maxT),rmslow(maxT), & rgylow(maxT),signlow(maxT),angnatlow(maxT), & qnatlow(maxT),rmsnatlow(maxT),rgynatlow(maxT), & angcorlow(maxT),qcorlow(maxT),rmscorlow(maxT),rgycorlow(maxT), & signcorlow(maxT),angcornatlow(maxT),qcornatlow(maxT), & rmscornatlow(maxT),rgycornatlow(maxT),signcornatlow(maxT), & angnormup(maxT),qup(maxT),rmsup(maxT),rgyup(maxT),signup(maxT), & angnatup(maxT),qnatup(maxT),rmsnatup(maxT),rgynatup(maxT), & angcorup(maxT),qcorup(maxT),rmscorup(maxT),rgycorup(maxT), & signcorup(maxT), & angcornatup(maxT),qcornatup(maxT),rmscornatup(maxT), & rgycornatup(maxT),signcornatup(MaxT) integer ilen,iroof external ilen,iroof double precision ebia(maxprot),rjunk character*10 liczba character*64 zeros / &'0000000000000000000000000000000000000000000000000000000000000000' & / logical lerr c print *,"Processor",me," calls read_protein_data" C Read seventh record: general data of proteins used for calibration call card_concat(controlcard,.true.) c write(2, *)controlcard call readi(controlcard,"NPROT",nprot,0) pdbref=index(controlcard,"PDBREF").gt.0 print_refstr=index(controlcard,"PRINT_REFSTR").gt.0 if (nprot.eq.0) then write(iout,*) "Error: at least one protein must be specified." call flush(iout) return1 endif protname(0)="all" do iprot=1,nprot read (inp,'(a)') protname(iprot) write (iout,*) write (iout,*) "Reading data of protein",iprot," named ", & protname(iprot) call card_concat(controlcard,.true.) call reada(controlcard,"ENECUT_MIN",enecut_min(iprot),15.0d0) call reada(controlcard,"ENECUT_MAX",enecut_max(iprot),100.0d0) call reada(controlcard,"ENECUT",enecut(iprot),enecut_min(iprot)) if (enecut(iprot).lt.enecut_min(iprot)) & enecut(iprot)=enecut_min(iprot) if (enecut_max(iprot).le.enecut_min(iprot)) & enecut_max(iprot)=2*enecut_min(iprot) write (iout,'(3(a,f9.1))') "ENECUT",enecut(iprot)," ENECUT_MIN", & enecut_min(iprot)," ENECUT_MAX",enecut_max(iprot) call readi(controlcard,"HEFAC",hefac(iprot),50) call readi(controlcard,"HTFAC",htfac(iprot),50) call readi(controlcard,"HELOW",hemax_low(iprot),20) call readi(controlcard,"HTLOW",htmax_low(iprot),20) write (iout,*) "iprot",iprot, & " hefac",hefac(iprot)," helow",hemax_low(iprot), & " htfac",htfac(iprot)," htlow",htmax_low(iprot) c 7/27/2013 AL Read maximum likelihood data call card_concat(controlcard,.true.) call readi(controlcard,"NBETA_L",nbeta(1,iprot),0) write (iout,*) "NBETA_L",nbeta(1,iprot) caonly(iprot)=index(controlcard,"CAONLY").gt.0 sconly(iprot)=index(controlcard,"SCONLY").gt.0 rmscomp(iprot)=index(controlcard,"RMSCOMP").gt.0 anglecomp(iprot)=index(controlcard,"ANGLECOMP").gt.0 sidecomp(iprot)=index(controlcard,"SIDECOMP").gt.0 call reada(controlcard,"SIGMA",sigma2(iprot),4.0d0) call reada(controlcard,"SIGMAANG",sigmaang2(iprot),4.0d0) call reada(controlcard,"SIGMASIDE",sigmaside2(iprot),4.0d0) write (iout,*) "RMSCOMP",rmscomp(iprot)," SIGMA",sigma2(iprot), & " CAONLY ",caonly(iprot)," SCONLY",sconly(iprot) write (iout,*) "ANGLECOMP",anglecomp(iprot), & " SIGMAANG",sigmaang2(iprot) write (iout,*) "SIDECOMP",sidecomp(iprot), & " SIGMASIDE",sigmaside2(iprot) do ib=1,nbeta(1,iprot) read(inp,*)betaT(ib,1,iprot),weilik(ib,iprot), & pdbfile(ib,iprot) write (iout,*) i,betaT(ib,1,iprot),weilik(ib,iprot), & pdbfile(ib,iprot) enddo c 7/27/2013 AL Read heat-capacity data call card_concat(controlcard,.true.) call readi(controlcard,"NBETA_CV",nbeta(2,iprot),0) call reada(controlcard,"WCV",wcv(iprot),1.0d0) call reada(controlcard,"BASE",heatbase(iprot),0.0d0) write (iout,*) "NBETA_CV",nbeta(2,iprot)," WCV",wcv(iprot) do ib=1,nbeta(2,iprot) read(inp,*) betaT(ib,2,iprot),target_cv(ib,iprot), & weiCv(ib,iprot) weiCv(ib,iprot)=weiCv(ib,iprot)*wcv(iprot) write (iout,*) betaT(ib,2,iprot),target_cv(ib,iprot), & weiCv(ib,iprot) enddo write (iout,*) "Experimental heat capacity curve" do ib=1,nbeta(2,iprot) write (iout,'(f6.2,2f10.5)') betaT(ib,2,iprot), & target_cv(ib,iprot),weiCv(ib,iprot) enddo write (iout,'(a,f10.5)') "Baseline",heatbase(iprot) call flush(iout) c Conformation-dependent averages call card_concat(controlcard,.true.) call readi(controlcard,"NATLIKE",natlike(iprot),0) do i=1,natlike(iprot) call card_concat(controlcard,.true.) call reada(controlcard,"WNAT",wnat(i,iprot),1.0d0) call readi(controlcard,"NUMNAT",numnat(i,iprot),1) call readi(controlcard,"NATDIM",natdim(i,iprot),1) do ib=1,nbeta(i+2,iprot) read(inp,*)betaT(ib,i+2,iprot),(weinat(k,ib,i,iprot), & nuexp(k,ib,i,iprot),k=1,natdim(i,iprot)) enddo enddo do i=1,natlike(iprot)+2 do j=1,nbeta(i,iprot) betaT(j,i,iprot)=1.0d0/(Rgas*betaT(j,i,iprot)) write (2,*) "R i",i," j",j," beta",betaT(j,i,iprot) enddo enddo do iparm=1,nparmset C Read names of files with the data for protein IPROT call card_concat(controlcard,.false.) if (iparm.eq.myparm) then call split_string(controlcard,protfiles(1,iprot), & maxfile_prot,nfile_prot(iprot)) #ifdef DEBUG write(iout,*)"iprot",iprot," nfile",nfile_prot(iprot) write(iout,*) & (protfiles(i,iprot),i=1,nfile_prot(iprot)) #endif endif enddo ! iparm c Read molecular information of the protein call molread_zs(iprot) c 3/31/04 AL Read the reference structures of protein iprot c print *,"Calling read_ref_structure" call read_ref_structure(iprot,*11) #ifdef DEBUG write (iout,*) "Protein",iprot," left READ_REF_STRUCTURE" #endif enddo ! iprot return 11 return1 end c------------------------------------------------------------------------------- subroutine read_database(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" #endif include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.PROTNAME" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.ALLPROT" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.VAR" include "COMMON.SBRIDGE" include "COMMON.GEO" include "COMMON.COMPAR" include "COMMON.OPTIM" character*128 nazwa character*16000 controlcard character*16000 all_protfiles character*4 liczba,liczba1 integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,if,ib,ibatch, & nn,nn1,inan integer ixdrf,iret,itmp integer nrec,nlines,iscor double precision energ,t_acq,tcpu integer ilen,iroof external ilen,iroof double precision rjunk integer ntot_all(0:maxprot,0:maxprocs-1) logical lerr double precision energia(0:max_ene),etot real*4 csingle(3,maxres2),reini,refree,rmsdev,prec integer Previous,Next c print *,"Processor",me," calls read_protein_data" #ifdef MPI if (me.eq.master) then Previous=MPI_PROC_NULL else Previous=me-1 endif if (me.eq.nprocs-1) then Next=MPI_PROC_NULL else Next=me+1 endif #endif c Set the scratchfile names write (liczba,'(bz,i4.4)') me write (liczba1,'(bz,i4.4)') myparm do iprot=1,nprot c 1/27/05 AL Change stored coordinates to single precision and don't store c energy components in the binary databases. lenrec(iprot)=12*(nres_zs(iprot)+nct_zs(iprot)-nnt_zs(iprot)+1) & +4*(2*nss_zs(1,iprot)+1)+8*natlike(iprot)*maxdimnat+16 c 4/13/04 AL Add space for similarity measure lenrec_ene(iprot) = (2*nntyp+3*n_ene+2)*8 & +8*natlike(iprot)*maxdimnat #ifdef DEBUG call flush(iout) write (iout,*) "Protein i"," lenrec",lenrec(iprot) write (iout,*) "lenrec_ene",lenrec_ene(iprot) call flush(iout) #endif bprotfiles(iprot) = scratchdir(:ilen(scratchdir))// & "/"//protname(iprot)(:ilen(protname(iprot)))// & "_par"//liczba1//"_"//liczba//".xbin" benefiles(iprot) = scratchdir(:ilen(scratchdir))// & "/"//protname(iprot)(:ilen(protname(iprot)))// & "_par"//liczba1//"_"//liczba//".enebin" c write (iout,*) "scratchfile ", c & bprotfiles(iprot)(:ilen(bprotfiles(iprot))) enddo call flush(iout) do i=0,nprot ntot(i)=0 enddo do iprot=1,nprot jj_old=1 call restore_molinfo(iprot) open (ientout,file=bprotfiles(iprot),status="unknown", & form="unformatted",access="direct",recl=lenrec(iprot)) c Change AL 12/30/2017 if (.not.mod_other_params) & open (istat,file=benefiles(iprot),status="unknown", & form="unformatted",access="direct",recl=lenrec_ene(iprot)) c Read conformations from binary DA files (one per batch) and write them to c a binary DA scratchfile. jj=0 jjj=0 write (liczba,'(bz,i4.4)') me #ifdef MPI IF (ME.EQ.MASTER) THEN c Only the master reads the database; it'll send it to the other procs c through a ring. #endif t_acq = tcpu() icount=0 itmp=0 do if=1,nfile_prot(iprot) nazwa=protfiles(if,iprot) & (:ilen(protfiles(if,iprot)))//".cx" #ifdef DEBUG write (iout,*) "Opening file ",nazwa(:ilen(nazwa)) #endif #if (defined(AIX) && !defined(JUBL)) call xdrfopen_(ixdrf,nazwa, "r", iret) #else call xdrfopen(ixdrf,nazwa, "r", iret) #endif if (iret.eq.0) goto 1111 do i=1,maxstr prec=10000.0 #if (defined(AIX) && !defined(JUBL)) call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret) if (iret.eq.0) goto 102 call xdrfint_(ixdrf, nss, iret) if (iret.eq.0) goto 102 do j=1,nss call xdrfint_(ixdrf, ihpb(j), iret) if (iret.eq.0) goto 102 call xdrfint_(ixdrf, jhpb(j), iret) if (iret.eq.0) goto 102 enddo call xdrffloat_(ixdrf,reini,iret) if (iret.eq.0) goto 102 call xdrffloat_(ixdrf,refree,iret) if (iret.eq.0) goto 102 call xdrfint_(ixdrf,natlik,iret) if (iret.eq.0) goto 102 do j=1,natlik call xdrfint(ixdrf,natliktemp(j),iret) if (iret.eq.0) goto 102 do k=1,natliktemp(j) call xdrffloat(ixdrf,nutemp(k,j),iret) if (iret.eq.0) goto 102 enddo enddo #else c write (0,*) "me",me," iprot",iprot," i",i call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret) if (iret.eq.0) goto 102 call xdrfint(ixdrf, nss, iret) if (iret.eq.0) goto 102 do k=1,nss call xdrfint(ixdrf, ihpb(k), iret) if (iret.eq.0) goto 102 call xdrfint(ixdrf, jhpb(k), iret) if (iret.eq.0) goto 102 enddo call xdrffloat(ixdrf,reini,iret) if (iret.eq.0) goto 102 call xdrffloat(ixdrf,refree,iret) if (iret.eq.0) goto 102 #ifdef FINAL call xdrfint(ixdrf,natlik,iret) if (iret.eq.0) goto 102 do j=1,natlik call xdrfint(ixdrf,natliktemp(j),iret) if (iret.eq.0) goto 102 do k=1,natliktemp(j) call xdrffloat(ixdrf,nutemp(k,j),iret) if (iret.eq.0) goto 102 enddo enddo #else natlik=0 call xdrffloat(ixdrf,rmsdev,iret) if (iret.eq.0) goto 102 c write (2,*) "rmsdev",rmsdev call xdrfint(ixdrf,iscor,iret) if (iret.eq.0) goto 102 c write (2,*) "iscor",iscor #endif #endif eini(jj+1,iprot)=reini entfac(jj+1,iprot)=refree do k=1,nres do l=1,3 c(l,k)=csingle(l,k) enddo enddo do k=nnt,nct do l=1,3 c(l,nres+k)=csingle(l,nres+k-nnt+1) enddo enddo #ifdef DEBUG write (iout,'(5hREAD ,i5,2f15.4)') & jj+1,eini(jj+1,iprot),entfac(jj+1,iprot) write (iout,*) "natlik",natlik do l=1,natlik write (iout,*) "natliketemp",natliktemp(l) write(iout,*) (nutemp(k,l),k=1,natliktemp(l)) enddo write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) call flush(iout) #endif call add_new_cconf(jjj,jj,jj_old,icount,iprot, & Next) enddo 102 continue write (iout,*) "Protein ",protname(iprot), & i-1," conformations read from DA file ", & nazwa(:ilen(nazwa)) write (iout,*) jj," conformations read so far" #if (defined(AIX) && !defined(JUBL)) call xdrfclose_(ixdrf,nazwa,iret) #else call xdrfclose(ixdrf,nazwa,iret) #endif c print *,"file ",nazwa(:ilen(nazwa))," closed" enddo #ifdef MPI #ifdef DEBUG write (iout,*) "jj_old",jj_old," jj",jj #endif call write_and_send_cconf(icount,jj_old,jj,iprot,Next) if (icount.gt.0) call MPI_Send(0,1,MPI_INTEGER,Next,570, & WHAM_COMM,IERROR) jj_old=jj+1 #else call write_and_send_cconf(icount,jj_old,jj,iprot,Next) #endif t_acq = tcpu() - t_acq write (iout,*) "Processor",me," protein",iprot, & " batch",ibatch," time for conformation read/send",t_acq #ifdef MPI ELSE c A worker gets the confs from the master and sends them to its neighbor t_acq = tcpu() call receive_and_pass_cconf(icount,jj_old,jj,iprot, & Previous,Next) t_acq = tcpu() - t_acq write (iout,*) "Processor",me," protein",iprot, & " batch",ibatch, & " time for conformation receive/send",t_acq ENDIF #endif ntot(iprot)=jj write (iout,*) "Protein", & protname(iprot)(:ilen(protname(iprot))),", ",ntot(iprot), & " conformatons read ",ntot(iprot)," conformations written to ", & bprotfiles(iprot)(:ilen(bprotfiles(iprot))) ntot(0) = ntot(0)+ntot(iprot) close(ientout) close(istat) enddo ! iprot write(iout,*)"A total of",ntot(0)," conformations read." #ifdef MPI c Check if everyone has the same number of conformations call MPI_Allgather(ntot(0),maxprot+1,MPI_INTEGER, & ntot_all(0,0),maxprot+1,MPI_INTEGER,MPI_Comm_World,IERROR) lerr=.false. do i=0,nprocs-1 if (i.ne.me) then do j=1,nprot if (ntot(j).ne.ntot_all(j,i)) then write (iout,*) "Number of conformations at processor",i, & " for protein",j," differs from that at processor",me, & ntot(j),ntot_all(j,i) lerr = .true. endif enddo endif enddo #ifdef WEICLUST c----------- Temporary; reading probs from external file open (88,file='1LE1_wham_last_2.rms',status='old') do i=1,ntot(1) read (88,*) ii,weirms(i,1) enddo do i=1,ntot(1) write (iout,*) "i",i," weirms",weirms(i,1) enddo close(88) call MPI_Bcast(weirms(1,1), ntot(1), MPI_Double_Precision, & Master, MPI_COMM_WORLD, IERROR) c----------- end temportary stuff #endif if (lerr) then write (iout,*) write (iout,*) "Number of conformations read by processors" write (iout,'(10x,7a10)') (protname(i),i=0,nprot) write (iout,*) do i=0,nprocs-1 write (iout,'(8i10)') i,(ntot_all(j,i),j=0,nprot) enddo write (iout,*) "Calculation terminated." call flush(iout) return1 endif return #endif 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa)) call flush(iout) #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) #endif return1 end c------------------------------------------------------------------------------ subroutine add_new_cconf(jjj,jj,jj_old,icount,iprot,Next) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.LOCAL" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.PROTNAME" include "COMMON.VMCPAR" include "COMMON.WEIGHTS" include "COMMON.NAMES" include "COMMON.ALLPROT" include "COMMON.WEIGHTDER" include "COMMON.VAR" include "COMMON.SBRIDGE" include "COMMON.GEO" include "COMMON.COMPAR" #ifdef ISNAN logical isnan #endif integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,ib,ibatch, & nn,nn1,inan,Next,itj double precision etot,energia(0:max_ene) jjj=jjj+1 c if (entfac(jj+1,iprot).gt.-10.0d0 c & .or. entfac(jj+1,iprot).lt.-150.0d0) then c write (iout,*) "Entropy factor out of range for conformation", c & jjj,entfac(jj+1,iprot),", conformation skipped." c return c endif call int_from_cart1(.false.) do j=nnt+1,nct if (itype(j-1).eq.ntyp1 .or. itype(j).eq.ntyp1) cycle if (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.5d0) then write (iout,*) "nnt",nnt," nct",nct write (iout,*) "Bad CA-CA bond length",j," ",vbld(j) write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j) write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) & "This conformation WILL NOT be added to the database." return endif enddo do j=nnt,nct itj=itype(j) if (itype(j).ne.10 .and. itype(j).ne.ntyp1 & .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then write (iout,*) "nnt",nnt," nct",nct write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j) write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) & "This conformation WILL NOT be added to the database." return endif enddo do j=3,nres if (itype(j).eq.ntyp1 .or. itype(j-1).eq.ntyp1 & .or. itype(j-1).eq.ntyp1) cycle if (theta(j).le.0.0d0) then write (iout,*) & "Zero theta angle(s) in conformation",ii write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) & "This conformation WILL NOT be added to the database." return endif if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad enddo if (.not. init_ene) then etot=0.0d0 do j=1,n_ene etot=etot+ww(j)*enetb(icount+1,j,iprot) enddo inan=0 c detecting NaNQ #ifdef ISNAN #ifdef AIX if (isnan(etot).ne.0) inan=1 #else if (isnan(etot)) inan=1 #endif #else #ifdef WINPGI idumm=proc_proc(etot,inan) #else call proc_proc(etot,inan) #endif #endif else inan=0 endif if (inan.eq.1) then write (iout,*) "NaNs detected in some of the energy", & " components for protein",iprot," batch",ibatch, & " conformation",jjj write (iout,*) "etot",etot write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres) write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct) write (iout,*) "The internal geometry is:" write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) "The components of the energy are:" energia(0)=etot do k=1,n_ene energia(k)=enetb(jj+1,k,iprot) enddo call enerprint(energia(0)) write (iout,*) & "This conformation WILL NOT be added to the database." call flush(iout) else jj=jj+1 #ifdef DEBUG write (iout,'(e15.5,16i5)') entfac(icount+1,iprot), & iscore(icount+1,0,iprot),ibatch write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres) write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct) write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct) write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss) write (iout,'(5e15.5)') (enetb(icount+1,j,iprot),j=1,n_ene) c write (iout,'(2e15.5)') (eneps(1,j,icount+1,iprot), c & eneps(2,j,icount+1,iprot),j=1,nntyp) #endif c write (iout,*) "First nu in readrtms" call flush(iout) icount=icount+1 list_conf(jj,iprot)=jj call store_cconf_from_file(jj,icount,iprot) do j=1,natlike(iprot) do k=1,natlik if (k.eq.numnat(j,iprot)) then do l=1,natdim(j,iprot) nu(l,k,j,iprot)=nutemp(l,k) enddo exit endif enddo enddo c write (iout,*) "End First nu in readrtms" call flush(iout) if (icount.eq.maxstr_proc) then #ifdef DEBUG write (iout,* ) "jj_old",jj_old," jj",jj write (iout,*) "Calling write_and_send_cconf" call flush(iout) #endif call write_and_send_cconf(icount,jj_old,jj,iprot,Next) jj_old=jj+1 #ifdef DEBUG write (iout,*) "Exited write_and_send_cconf" call flush(iout) #endif icount=0 endif endif return end c------------------------------------------------------------------------------ subroutine store_cconf_from_file(jj,icount,iprot) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CHAIN" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.ALLPROT" include "COMMON.VAR" integer i,j,jj,icount,ibatch,iprot c Store the conformation that has been read in do i=1,2*nres do j=1,3 c_zs(j,i,icount,iprot)=c(j,i) enddo enddo nss_zs(icount,iprot)=nss do i=1,nss ihpb_zs(i,icount,iprot)=ihpb(i) jhpb_zs(i,icount,iprot)=jhpb(i) enddo return end c------------------------------------------------------------------------------ subroutine write_and_send_cconf(icount,jj_old,jj,iprot,Next) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR include "COMMON.MPI" #endif include "COMMON.WEIGHTS" include "COMMON.CHAIN" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.VAR" include "COMMON.ALLPROT" include "COMMON.ENERGIES" include "COMMON.WEIGHTDER" include "COMMON.OPTIM" include "COMMON.COMPAR" integer icount,jj_old,jj,Next,iprot c Write the structures to a scratch file #ifdef MPI c Master sends the portion of conformations that have been read in to the neighbor #ifdef DEBUG write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF" call flush(iout) #endif call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,IERROR) #ifdef DEBUG write (iout,*) "Processor",me," Next",next," sent icount=",icount write (iout,*) "Processor",me," jj_old",jj_old," jj",jj call flush(iout) #endif if (icount.gt.0) then call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER, & Next,571,WHAM_COMM,IERROR) call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER, & Next,572,WHAM_COMM,IERROR) if (.not. init_ene) call MPI_Send(enetb(jj_old,1,iprot), & maxstr*n_ene, & MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,IERROR) call MPI_Send(nu(1,1,jj_old,iprot), & maxdimnat*natlike(iprot)*icount, & MPI_DOUBLE_PRECISION, & Next,577,WHAM_COMM,IERROR) call MPI_Send(eini(jj_old,iprot),icount, & MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR) call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION, & Next,579,WHAM_COMM,IERROR) call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2, & MPI_REAL,Next,580,WHAM_COMM,IERROR) if (.not. init_ene) call MPI_Send(eneps(1,1,1,iprot), & 2*icount*nntyp, & MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR) endif #endif call flush(iout) call dawrite_ccoords(iprot,jj_old,jj,ientout) c Change AL 20/12/2017 if (.not. mod_other_params) &call dawrite_ene(iprot,jj_old,jj,istat) return end c------------------------------------------------------------------------------ #ifdef MPI subroutine receive_and_pass_cconf(icount,jj_old,jj,iprot,Previous, & Next) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "mpif.h" integer IERROR,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" include "COMMON.CHAIN" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.ALLPROT" include "COMMON.ENERGIES" include "COMMON.VAR" include "COMMON.GEO" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.COMPAR" include "COMMON.OPTIM" integer i,j,k,icount,jj_old,jj,iprot,Previous,Next icount=1 #ifdef DEBUG write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF" call flush(iout) #endif do while (icount.gt.0) c call MPI_Probe(Previous,570,WHAM_COMM,STATUS,IERROR) call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,WHAM_COMM, & STATUS,IERROR) #ifdef DEBUG write (iout,*)"Processor",me," previous",previous," icount",icount call flush(iout) #endif call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM, & IERROR) #ifdef DEBUG write (iout,*) "Processor",me," icount",icount call flush(iout) #endif if (icount.eq.0) return call MPI_Recv(nss_zs(1,iprot),icount,MPI_INTEGER, & Previous,571,WHAM_COMM,STATUS,IERROR) call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER, & Next,571,WHAM_COMM,IERROR) call MPI_Recv(ihpb_zs(1,1,iprot),icount,MPI_INTEGER, & Previous,572,WHAM_COMM,STATUS,IERROR) call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER, & Next,572,WHAM_COMM,IERROR) if (.not. init_ene) then call MPI_Recv(enetb(jj_old,1,iprot),maxstr*n_ene, & MPI_DOUBLE_PRECISION,Previous,576,WHAM_COMM,STATUS,IERROR) call MPI_Send(enetb(jj_old,1,iprot),maxstr*n_ene, & MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,STATUS,IERROR) endif call MPI_Recv(nu(1,1,jj_old,iprot), & maxdimnat*natlike(iprot)*icount, & MPI_DOUBLE_PRECISION, & Previous,577,WHAM_COMM,STATUS,IERROR) call MPI_Send(nu(1,1,jj_old,iprot), & maxdimnat*natlike(iprot)*icount, & MPI_DOUBLE_PRECISION, & Next,577,WHAM_COMM,IERROR) call MPI_Recv(eini(jj_old,iprot),icount, & MPI_DOUBLE_PRECISION,Previous,578,WHAM_COMM,STATUS,IERROR) call MPI_Send(eini(jj_old,iprot),icount, & MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR) call MPI_Recv(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION, & Previous,579,WHAM_COMM,STATUS,IERROR) call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION, & Next,579,WHAM_COMM,IERROR) call MPI_Recv(c_zs(1,1,1,iprot),3*icount*maxres2, & MPI_REAL,Previous,580,WHAM_COMM,STATUS,IERROR) call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2, & MPI_REAL,Next,580,WHAM_COMM,IERROR) if (.not. init_ene) then call MPI_Recv(eneps(1,1,1,iprot),2*icount*nntyp, & MPI_DOUBLE_PRECISION,Previous,582,WHAM_COMM,STATUS,IERROR) call MPI_Send(eneps(1,1,1,iprot),2*icount*nntyp, & MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR) endif jj=jj_old+icount-1 do i=jj_old,jj list_conf(i,iprot)=i enddo call dawrite_ccoords(iprot,jj_old,jj,ientout) c Change AL 12/20/2017 if (.not. mod_other_params) &call dawrite_ene(iprot,jj_old,jj,istat) jj_old=jj+1 #ifdef DEBUG write (iout,*) "Protein",iprot write (iout,*) "Processor",me," received",icount," conformations" do i=1,icount write (iout,'(8f10.4)') ((c_zs(l,k,i,iprot),l=1,3,k=1,nres) write (iout,'(8f10.4)')((c_zs(l,k,i+nres,iprot),l=1,3,k=nnt,nct) write (iout,'(16i5)') nss_zs(i,iprot),(ihpb_zs(k,i,iprot), & jhpb_zs(k,i,iprot),k=1,nss_zs(i,iprot)) write (iout,'(5e15.5)') (enetb(i,j,iprot),j=1,n_ene) write (iout,'(2e15.5)') (eneps(1,j,i,iprot), & eneps(2,j,i,iprot),j=1,nntyp) write (iout,'(e15.5,16i5)') entfac(i),iscore(i,0,iprot), & kbatch(i,iprot) enddo #endif enddo return end #endif c------------------------------------------------------------------------------ subroutine read_thermal implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.CHAIN" include "COMMON.SBRIDGE" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.CLASSES" include "COMMON.VAR" include "COMMON.THERMAL" character*800 controlcard call card_concat(controlcard,.true.) call readi(controlcard,"NGRIDT",NGridT,200) call reada(controlcard,"DELTAT",deltaT,5.0d0) call reada(controlcard,"T0",GridT0,2.0d2) write (iout,*) "Parameters of thermal curves" write (iout,*) "NGRIDT",NGridT," DELTAT",deltaT," T0",GridT0 return end c------------------------------------------------------------------------------ subroutine daread_ene(iprot,istart_conf,iend_conf) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" #endif include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.ALLPROT" include "COMMON.WEIGHTDER" include "COMMON.COMPAR" include "COMMON.VMCPAR" integer iprot,istart_conf,iend_conf integer i,ii,iii,j,k #ifdef DEBUG write (iout,*) "Calling DAREAD_ENE" #endif c write (iout,*) "Reading: n_ene",n_ene c call flush(iout) c write (iout,*) "DAREAD_ENE",istart_conf,iend_conf c c Read conformations off a DA scratchfile. c do ii=istart_conf,iend_conf iii=list_conf(ii,iprot) i = ii - istart_conf + 1 read(ientin,rec=iii) (enetb(i,j,iprot),j=1,n_ene), & (enetb_orig(i,j,iprot),j=1,n_ene), & (enetb_oorig(i,j,iprot),j=1,n_ene), & eini(ii,iprot),entfac(ii,iprot), & (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp), & ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot)) #ifdef DEBUG write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot), & entfac(ii,iprot), write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene) write (iout,'(18(1pe12.4))') & ((nu(k,j,i,iprot)k=1,maxdimnat),j=1,natlike(iprot)) call flush(iout) #endif enddo return end c------------------------------------------------------------------------------ subroutine dawrite_ene(iprot,istart_conf,iend_conf,unit_out) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" #endif include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.ALLPROT" include "COMMON.WEIGHTDER" include "COMMON.VMCPAR" include "COMMON.COMPAR" integer iprot,istart_conf,iend_conf,unit_out integer i,ii,iii,j,k c write (iout,*) "Writing: n_ene",n_ene c call flush(iout) c write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf c c Write conformations to a DA scratchfile. c do ii=istart_conf,iend_conf iii=list_conf(ii,iprot) i = ii - istart_conf + 1 write(unit_out,rec=iii) (enetb(i,j,iprot),j=1,n_ene), & (enetb_orig(i,j,iprot),j=1,n_ene), & (enetb_oorig(i,j,iprot),j=1,n_ene), & eini(ii,iprot),entfac(ii,iprot), & (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp), & ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot)) #ifdef DEBUG write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot), & entfac(ii,iprot) write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene) write (iout,'(18(1pe12.4))') & ((nu(kj,i,iprot),k=1,maxdimnat),j=1,natlike(iprot)) call flush(iout) #endif enddo return end c------------------------------------------------------------------------------ subroutine daread_ccoords(iprot,istart_conf,iend_conf) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" #endif include "COMMON.CHAIN" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.ALLPROT" include "COMMON.INTERACT" include "COMMON.VAR" include "COMMON.SBRIDGE" include "COMMON.GEO" include "COMMON.COMPAR" include "COMMON.VMCPAR" include "COMMON.WEIGHTDER" integer iprot,istart_conf,iend_conf integer i,j,k,ij,ii,iii integer len real*4 csingle(3,maxres2) character*16 form,acc character*32 nam c c Read conformations off a DA scratchfile. c #ifdef DEBUG write (iout,*) "DAREAD_COORDS" write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf write (iout,*) "lenrec",lenrec(iprot) inquire(unit=ientin,name=nam,recl=len,form=form,access=acc) write (iout,*) "len=",len," form=",form," acc=",acc write (iout,*) "nam=",nam call flush(iout) #endif do ii=istart_conf,iend_conf iii=list_conf(ii,iprot) ij = ii - istart_conf + 1 #ifdef DEBUG write (iout,*) "Reading binary file, record",iii," ii",ii call flush(iout) #endif read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres), & ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres), & nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot), & entfac(ii,iprot), & ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot)) do i=1,2*nres do j=1,3 c(j,i)=csingle(j,i) enddo enddo #ifdef DEBUG write (iout,*) "iprot",iprot," ii",ii write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres) write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres) write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot) write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss) call flush(iout) #endif call store_ccoords(iprot,ii-istart_conf+1) enddo return end c------------------------------------------------------------------------------ subroutine dawrite_ccoords(iprot,istart_conf,iend_conf,unit_out) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" include "COMMON.MPI" #endif include "COMMON.CHAIN" include "COMMON.INTERACT" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.ALLPROT" include "COMMON.VAR" include "COMMON.SBRIDGE" include "COMMON.GEO" include "COMMON.COMPAR" include "COMMON.WEIGHTDER" include "COMMON.VMCPAR" real*4 csingle(3,maxres2) integer iprot,istart_conf,iend_conf integer i,j,k,ii,ij,iii,unit_out integer len character*16 form,acc character*32 nam c c Write conformations to a DA scratchfile. c #ifdef DEBUG write (iout,*) "DAWRITE_COORDS" write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf write (iout,*) "lenrec",lenrec(iprot) inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc) write (iout,*) "len=",len," form=",form," acc=",acc write (iout,*) "nam=",nam call flush(iout) #endif do ii=istart_conf,iend_conf iii=list_conf(ii,iprot) ij = ii - istart_conf + 1 call restore_ccoords(iprot,ii-istart_conf+1) #ifdef DEBUG write (iout,*) "Writing binary file, record",iii," ii",ii call flush(iout) #endif do i=1,2*nres do j=1,3 csingle(j,i)=c(j,i) enddo enddo write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres), & ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres), & nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot), & entfac(ii,iprot), & ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot)) #ifdef DEBUG write (iout,*) "kbatch",kbatch(ii,iprot)," iscore", & iscore(ii,0,iprot) write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres) write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres) write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot) write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss) call flush(iout) #endif enddo return end c------------------------------------------------------------------------------ subroutine store_ccoords(iprot,ii) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.VAR" include "COMMON.CHAIN" include "COMMON.ALLPROT" include "COMMON.IOUNITS" include "COMMON.GEO" include "COMMON.SBRIDGE" double precision thetnorm integer iprot,i,j,ii do i=1,nres_zs(iprot) do j=1,3 c_zs(j,i,ii,iprot)=c(j,i) enddo enddo do i=nnt_zs(iprot),nct_zs(iprot) do j=1,3 c_zs(j,i+nres,ii,iprot)=c(j,i+nres) enddo enddo c 5/7/02 AL: store sbridge info nss_zs(ii,iprot)=nss do i=1,nss ihpb_zs(i,ii,iprot)=ihpb(i) jhpb_zs(i,ii,iprot)=jhpb(i) enddo return end c------------------------------------------------------------------------------ subroutine restore_ccoords(iprot,ii) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.INTERACT" include "COMMON.VAR" include "COMMON.ALLPROT" include "COMMON.SBRIDGE" include "COMMON.CHAIN" include "COMMON.IOUNITS" integer iprot,i,j,ii do i=1,nres_zs(iprot) do j=1,3 c(j,i)=c_zs(j,i,ii,iprot) enddo enddo do i=nnt_zs(iprot),nct_zs(iprot) do j=1,3 c(j,i+nres)=c_zs(j,i+nres,ii,iprot) enddo enddo c 5/7/02 AL: restore sbridge info nss=nss_zs(ii,iprot) do i=1,nss ihpb(i)=ihpb_zs(i,ii,iprot)+nres jhpb(i)=jhpb_zs(i,ii,iprot)+nres c forcon(i)=fbr c dhpb(i)=dbr enddo #ifdef DEBUG write (iout,*) "restore_ccoords",ii write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres) write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres) write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss) call flush(iout) #endif return end c------------------------------------------------------------------------------ subroutine card_concat(card,to_upper) implicit none include 'DIMENSIONS.ZSCOPT' include "COMMON.IOUNITS" character*(*) card character*80 karta,ucase logical to_upper integer ilen external ilen read (inp,'(a)') karta if (to_upper) karta=ucase(karta) card=' ' do while (karta(80:80).eq.'&') card=card(:ilen(card)+1)//karta(:79) read (inp,'(a)') karta if (to_upper) karta=ucase(karta) enddo card=card(:ilen(card)+1)//karta return end c------------------------------------------------------------------------------ subroutine readi(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch integer wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end c---------------------------------------------------------------------------- subroutine reada(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch character*80 aux double precision wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end c---------------------------------------------------------------------------- subroutine multreadi(rekord,lancuch,tablica,dim,default) implicit none integer dim,i integer tablica(dim),default character*(*) rekord,lancuch character*80 aux integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) return iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end c---------------------------------------------------------------------------- subroutine multreada(rekord,lancuch,tablica,dim,default) implicit none integer dim,i double precision tablica(dim),default character*(*) rekord,lancuch character*80 aux integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) return iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end c---------------------------------------------------------------------------- subroutine reads(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch,wartosc,default character*80 aux integer ilen,lenlan,lenrec,iread,ireade external ilen logical iblnk external iblnk lenlan=ilen(lancuch) lenrec=ilen(rekord) iread=index(rekord,lancuch(:lenlan)//"=") c print *,"rekord",rekord," lancuch",lancuch c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec if (iread.eq.0) then wartosc=default return endif iread=iread+lenlan+1 c print *,"iread",iread c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) do while (iread.le.lenrec .and. iblnk(rekord(iread:iread))) iread=iread+1 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) enddo c print *,"iread",iread if (iread.gt.lenrec) then wartosc=default return endif ireade=iread+1 c print *,"ireade",ireade do while (ireade.lt.lenrec .and. & .not.iblnk(rekord(ireade:ireade))) ireade=ireade+1 enddo wartosc=rekord(iread:ireade) return end c---------------------------------------------------------------------------- subroutine multreads(rekord,lancuch,tablica,dim,default) implicit none integer dim,i character*(*) rekord,lancuch,tablica(dim),default character*80 aux integer ilen,lenlan,lenrec,iread,ireade external ilen logical iblnk external iblnk do i=1,dim tablica(i)=default enddo lenlan=ilen(lancuch) lenrec=ilen(rekord) iread=index(rekord,lancuch(:lenlan)//"=") c print *,"rekord",rekord," lancuch",lancuch c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec if (iread.eq.0) return iread=iread+lenlan+1 do i=1,dim c print *,"iread",iread c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) do while (iread.le.lenrec .and. iblnk(rekord(iread:iread))) iread=iread+1 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) enddo c print *,"iread",iread if (iread.gt.lenrec) return ireade=iread+1 c print *,"ireade",ireade do while (ireade.lt.lenrec .and. & .not.iblnk(rekord(ireade:ireade))) ireade=ireade+1 enddo tablica(i)=rekord(iread:ireade) iread=ireade+1 enddo end c---------------------------------------------------------------------------- subroutine split_string(rekord,tablica,dim,nsub) implicit none integer dim,nsub,i,ii,ll,kk character*(*) tablica(dim) character*(*) rekord integer ilen external ilen do i=1,dim tablica(i)=" " enddo ii=1 ll = ilen(rekord) nsub=0 do i=1,dim C Find the start of term name kk = 0 do while (ii.le.ll .and. rekord(ii:ii).eq." ") ii = ii+1 enddo C Parse the name into TABLICA(i) until blank found do while (ii.le.ll .and. rekord(ii:ii).ne." ") kk = kk+1 tablica(i)(kk:kk)=rekord(ii:ii) ii = ii+1 enddo if (kk.gt.0) nsub=nsub+1 if (ii.gt.ll) return enddo return end c------------------------------------------------------------------------------- integer function iroof(n,m) ii = n/m if (ii*m .lt. n) ii=ii+1 iroof = ii return end