module csa !----------------------------------------------------------------------------- use io_units use names use math use random use geometry_data, only: nres,rad2deg use geometry use energy use MPI_ use csa_data use compare use io_config implicit none !----------------------------------------------------------------------------- ! commom.bank ! common/varin/ ! real(kind=8),dimension(:,:,:,:),allocatable :: dihang_in !(mxang,maxres,mxch,mxio) integer,dimension(:),allocatable :: nss_in !(mxio) integer,dimension(:,:),allocatable :: iss_in,jss_in !(maxss,mxio) ! common/minvar/ real(kind=8),dimension(:,:,:,:),allocatable :: dihang !(mxang,maxres,mxch,mxio) real(kind=8),dimension(:),allocatable :: etot!,rmsn,pncn !(mxio) integer,dimension(:),allocatable :: nss_out !(mxio) integer,dimension(:,:),allocatable ::iss_out,jss_out !(maxss,mxio) ! common/bank/ ! real(kind=8),dimension(:,:,:,:),allocatable :: bvar !(mxang,maxres,mxch,mxio) ! real(kind=8),dimension(:),allocatable :: bene,rene,& ! brmsn,rrmsn,bpncn,rpncn !(mxio) integer,dimension(:),allocatable :: is,jbank !(mxio) real(kind=8) :: avedif,difmin,ebmin,ebmax,ebmaxt!,& ! dele,difcut,cutdif,rmscut,pnccut real(kind=8),dimension(:,:),allocatable :: dij !(mxio,mxio) ! common/bank_disulfid/ ! common/mvstat/ integer,dimension(:),allocatable :: movenx,movernx !(mxio) integer,dimension(:,:),allocatable :: nstatnx,nstatnx_tot !(0:mxmv,3) integer,dimension(:,:),allocatable :: indb !(mxio,9) integer,dimension(:,:),allocatable :: parent !(3,mxio) ! common/send2/ integer,dimension(:),allocatable :: isend2 !(mxio) integer,dimension(:,:),allocatable :: iff_in !(maxres,mxio2) integer,dimension(:,:,:,:),allocatable :: dihang_in2 !(mxang,maxres,mxch,mxio2) integer,dimension(:,:),allocatable :: idata !(5,mxio) !----------------------------------------------------------------------------- ! common.csa ! integer :: irestart,ndiff ! common/alphaa/ integer,dimension(:),allocatable :: ngroup !(mxgr) integer,dimension(:,:,:),allocatable :: igroup !(3,mxang,mxgr) integer :: ntotgr!,numch ! common/csa_input/ ! common/dih_control/ ! real(kind=8) :: rdih_bias logical :: ldih_bias !----------------------------------------------------------------------------- ! common.distfit ! COMMON /frag/ integer,dimension(:,:),allocatable :: bvar_frag !(mxio,6) integer,dimension(:,:),allocatable :: hvar_frag,lvar_frag,svar_frag !(mxio,3) integer,dimension(:,:),allocatable :: avar_frag !(mxio,5) !----------------------------------------------------------------------------- ! commom.hairpin ! common /spinka/ integer :: nharp_tot integer,dimension(:),allocatable :: nharp_seed,nharp_use !(max_seed) integer,dimension(:,:,:),allocatable :: iharp_seed !(4,maxres/3,max_seed) integer,dimension(:,:,:),allocatable :: iharp_use !(0:4,maxres/3,max_seed) !----------------------------------------------------------------------------- ! Maximum number of moves (n1-n8) integer,parameter :: mxmv=18 !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- ! bank.F !----------------------------------------------------------------------------- subroutine refresh_bank(ntrial) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.CONTROL' character :: chacc integer :: iaccn,ntrial real(kind=8) :: l_diff(mxio),denep integer :: i,j,n,m,i1,idmin real(kind=8) :: del_ene do i=0,mxmv do j=1,3 nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j) nstatnx(i,j)=0 enddo enddo ! loop over all newly obtained conformations do n=1,ntrial chacc=' ' iaccn=0 nstatnx(movernx(n),1)=nstatnx(movernx(n),1)+1 !ccccccccccccccccccccccccccccccccccccccccccc !jlee if(iref.ne.0) then if(rmsn(n).gt.rmscut.or.pncn(n).lt.pnccut) goto 100 endif !jlee if(etot(n).gt.ebmax) goto 100 ! Find the conformation closest to the conformation n in the bank difmin=9.d9 do m=1,nbank call get_diff12(dihang(1,1,1,n),bvar(1,1,1,m),l_diff(m)) if(l_diff(m).lt.difmin) then difmin=l_diff(m) idmin=m endif enddo if(difmin.lt.cutdif) then ! n is redundant to idmin if(etot(n).lt.bene(idmin)) then if(etot(n).lt.bene(idmin)-0.01d0) then ibank(idmin)=0 jbank(idmin)=0 endif denep=bene(idmin)-etot(n) call replace_bvar(idmin,n) !rc Update dij do i1=1,nbank if (i1.ne.idmin) then dij(i1,idmin)=l_diff(i1) dij(idmin,i1)=l_diff(i1) endif enddo chacc='c' iaccn=idmin nstatnx(movernx(n),2)=nstatnx(movernx(n),2)+1 if(idmin.eq.ibmax) call find_max endif else ! got new conformation del_ene=0.0d0 if(ebmax-ebmin.gt.del_ene) then denep=ebmax-etot(n) call replace_bvar(ibmax,n) !rc Update dij do i1=1,nbank if (i1.ne.ibmax) then dij(i1,ibmax)=l_diff(i1) dij(ibmax,i1)=l_diff(i1) endif enddo chacc='f' iaccn=ibmax nstatnx(movernx(n),3)=nstatnx(movernx(n),3)+1 ibank(ibmax)=0 jbank(ibmax)=0 call find_max else if(del_ene.lt.0.0001) then write (iout,*) 'ERROR in refresh_bank: ' write (iout,*) 'ebmax: ',ebmax write (iout,*) 'ebmin: ',ebmin write (iout,*) 'del_ene: ',del_ene !rc call mpi_abort(mpi_comm_world,ierror,ierrcode) endif !jp nbmax is never defined so condition below is always false ! if(nbank.lt.nbmax) then ! nbank=nbank+1 ! call replace_bvar(nbank,n) ! ibank(nbank)=0 ! jbank(nbank)=0 ! else call replace_bvar(ibmax,n) ibank(ibmax)=0 jbank(ibmax)=0 call find_max ! endif endif endif !ccccccccccccccccccccccccccccccccccccccccccc 100 continue if (iaccn.eq.0) then if (iref.eq.0) then write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5)') & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',& indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9) else write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,0pf4.1,a5,f3.0)') & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',& indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),& ' rms ',rmsn(n),' %NC ',pncn(n)*100 endif else if (iref.eq.0) then write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,1x,a1,i4,0pf8.1,0pf8.1)') & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',& indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),& chacc,iaccn,difmin,denep else write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4,3i5,a5,0pf4.1,a5,f3.0,1x,a1,i4,0pf8.1,0pf8.1)') & indb(n,2),' e ',indb(n,3),indb(n,1),' etot ',etot(n),' mv ',& indb(n,5),indb(n,4),indb(n,7),indb(n,8),indb(n,9),& ' rms ',rmsn(n),' %NC ',pncn(n)*100,& chacc,iaccn,difmin,denep endif endif enddo ! end of loop over all newly obtained conformations do i=0,mxmv if(nstatnx(i,1).ne.0) then if (i.le.9) then write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') & '## N',i,' total=',nstatnx(i,1),& ' close=',nstatnx(i,2),' far=',nstatnx(i,3),& ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1) else write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') & '##N',i,' total=',nstatnx(i,1),& ' close=',nstatnx(i,2),' far=',nstatnx(i,3),& ' %acc',(nstatnx(i,2)+nstatnx(i,3))*100.0/nstatnx(i,1) endif else if (i.le.9) then write(iout,'(a4,i1,a7,i4,a7,i4,a5,i4,a5,f5.1)') & '## N',i,' total=',nstatnx(i,1),& ' close=',nstatnx(i,2),' far=',nstatnx(i,3),& ' %acc',0.0 else write(iout,'(a3,i2,a7,i4,a7,i4,a5,i4,a5,f5.1)') & '##N',i,' total=',nstatnx(i,1),& ' close=',nstatnx(i,2),' far=',nstatnx(i,3),& ' %acc',0.0 endif endif enddo call flush(iout) !rc Update dij !rc moved up, saves some get_diff12 calls !rc !rc do i1=1,nbank-1 !rc do i2=i1+1,nbank !rc if(jbank(i1).eq.0.or.jbank(i2).eq.0) then !rc call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff) !rc dij(i1,i2)=diff !rc dij(i2,i1)=diff !rc endif !rc enddo !rc enddo do i=1,nbank jbank(i)=1 enddo return end subroutine refresh_bank !----------------------------------------------------------------------------- subroutine replace_bvar(iold,inew) use control_data, only: vdisulf use energy_data, only: ns,iss ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' integer :: iold,inew,ierror,ierrcode,i,j,k if (iold.gt.mxio .or. iold.lt.1 .or. inew.gt.mxio .or. inew.lt.1) & then write (iout,*) 'Dimension ERROR in REPLACE_BVAR: IOLD',iold,& ' INEW',inew call mpi_abort(mpi_comm_world,ierror,ierrcode) endif do k=1,numch do j=2,nres-1 do i=1,4 bvar(i,j,k,iold)=dihang(i,j,k,inew) enddo enddo enddo bene(iold)=etot(inew) brmsn(iold)=rmsn(inew) bpncn(iold)=pncn(inew) if(bene(iold).lt.ebmin) then ebmin=bene(iold) ibmin=iold endif if(vdisulf) then bvar_nss(iold)=nss_out(inew) !d write(iout,*) 'SS BANK',iold,bvar_nss(iold) do i=1,bvar_nss(iold) bvar_ss(1,i,iold)=iss_out(i,inew) bvar_ss(2,i,iold)=jss_out(i,inew) !d write(iout,*) 'SS',bvar_ss(1,i,iold)-nres, !d & bvar_ss(2,i,iold)-nres enddo bvar_ns(iold)=ns-2*bvar_nss(iold) !d write(iout,*) 'CYS #free ', bvar_ns(iold) k=0 do i=1,ns j=1 do while( iss(i).ne.iss_out(j,inew)-nres .and. & iss(i).ne.jss_out(j,inew)-nres .and. & j.le.nss_out(inew)) j=j+1 enddo if (j.gt.nss_out(inew)) then k=k+1 bvar_s(k,iold)=iss(i) endif enddo !d write(iout,*) 'CYS free',(bvar_s(k,iold),k=1,bvar_ns(iold)) endif return end subroutine replace_bvar !----------------------------------------------------------------------------- subroutine save_is(ind) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' integer :: ind,i,j,k,index,ierror,ierrcode index=nbank+ind ! print *, "nbank,ind,index,is(ind) ",nbank,ind,index,is(ind) if (index.gt.mxio .or. index.lt.1 .or. & is(ind).gt.mxio .or. is(ind).lt.1) then write (iout,*) 'Dimension ERROR in SAVE_IS: INDEX',index,& ' IND',ind,' IS',is(ind) call mpi_abort(mpi_comm_world,ierror,ierrcode) endif do k=1,numch do j=2,nres-1 do i=1,4 bvar(i,j,k,index)=bvar(i,j,k,is(ind)) enddo enddo enddo bene(index)=bene(is(ind)) ibank(is(ind))=1 return end subroutine save_is !----------------------------------------------------------------------------- subroutine select_is(n,ifar,idum) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer,dimension(mxio) :: itag real(kind=8),dimension(mxio) :: adiff integer :: n,ifar,idum,i,iusesv,imade iuse=0 do i=1,nbank if(ibank(i).eq.0) then iuse=iuse+1 itag(iuse)=i endif enddo iusesv=iuse if(iuse.eq.0) then icycle=icycle+1 do i=1,nbank if(ibank(i).eq.2) then ibank(i)=1 else ibank(i)=0 endif enddo imade=0 call get_is(idum,ifar,n,imade,0) !test3 call get_is_max(idum,ifar,n,imade,0) else if(iuse.eq.n) then do i=1,iuse is(i)=itag(i) call save_is(i) enddo else if(iuse.lt.n) then ! if(icycle.eq.0) then ! do i=1,n ! ind=mod(i-1,iuse)+1 ! is(i)=itag(ind) ! call save_is(i) ! enddo ! else ! endif do i=1,iuse is(i)=itag(i) call save_is(i) enddo imade=iuse ! call get_is_ran(idum,n,imade,1) call get_is(idum,ifar,n,imade,1) !test3 call get_is_max(idum,ifar,n,imade,1) ! if(iusesv.le.n/10) then if(iusesv.le.0) then icycle=icycle+1 do i=1,nbank ! if(ibank(i).eq.2) then ! ibank(i)=1 if(ibank(i).ge.2) then ibank(i)=ibank(i)-1 else ibank(i)=0 endif enddo endif else imade=0 call get_is(idum,ifar,n,imade,0) !test3 call get_is_max(idum,ifar,n,imade,0) endif iuse=iusesv return end subroutine select_is !----------------------------------------------------------------------------- subroutine get_is_ran(idum,n,imade,k) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! real(kind=4) :: ran1,ran2 integer,dimension(mxio) :: itag real(kind=8),dimension(mxio) :: adiff integer :: idum,n,imade,k,j,i,iran do j=imade+1,n iuse=0 do i=1,nbank if(ibank(i).eq.k) then iuse=iuse+1 itag(iuse)=i endif enddo iran=iuse* ran1(idum)+1 is(j)=itag(iran) call save_is(j) enddo return end subroutine get_is_ran !----------------------------------------------------------------------------- subroutine get_is(idum,ifar,n,imade,k) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! real(kind=4) :: ran1,ran2 integer,dimension(mxio) :: itag real(kind=8),dimension(mxio) :: adiff integer :: idum,ifar,n,imade,k,i,iran iuse=0 do i=1,nbank if(ibank(i).eq.k) then iuse=iuse+1 itag(iuse)=i endif enddo iran=iuse* ran1(idum)+1 imade=imade+1 is(imade)=itag(iran) call save_is(imade) do i=imade+1,ifar-1 if(icycle.eq.-1) then call select_iseed_max(i,k) else call select_iseed_min(i,k) !test4 call select_iseed_max(i,k) endif call save_is(i) enddo do i=ifar,n call select_iseed_far(i,k) call save_is(i) enddo return end subroutine get_is !----------------------------------------------------------------------------- subroutine select_iseed_max(imade1,ik) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer,dimension(mxio) :: itag real(kind=8),dimension(mxio) :: adiff integer :: imade1,ik,i,n,imade,m,itagi real(kind=8) :: difmax,diff,emax,benei,diffmn iuse=0 avedif=0.d0 difmax=0.d0 do n=1,nbank if(ibank(n).eq.ik) then iuse=iuse+1 diffmn=9.d190 do imade=1,imade1-1 ! m=nbank+imade ! call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff) m=is(imade) diff=dij(n,m) if(diff.lt.diffmn) diffmn=diff enddo if(diffmn.gt.difmax) difmax=diffmn adiff(iuse)=diffmn itag(iuse)=n avedif=avedif+diffmn endif enddo avedif=avedif/iuse ! avedif=(avedif+difmax)/2 emax=-9.d190 do i=1,iuse if(adiff(i).ge.avedif) then itagi=itag(i) benei=bene(itagi) if(benei.gt.emax) then emax=benei is(imade1)=itagi endif endif enddo if(ik.eq.0) iuse=iuse-1 return end subroutine select_iseed_max !----------------------------------------------------------------------------- subroutine select_iseed_min(imade1,ik) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer,dimension(mxio) :: itag real(kind=8),dimension(mxio) :: adiff integer :: imade1,ik,n,imade,m,i,itagi real(kind=8) :: difmax,diff,diffmn,emin,benei iuse=0 avedif=0.d0 difmax=0.d0 do n=1,nbank if(ibank(n).eq.ik) then iuse=iuse+1 diffmn=9.d190 do imade=1,imade1-1 ! m=nbank+imade ! call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff) m=is(imade) diff=dij(n,m) if(diff.lt.diffmn) diffmn=diff enddo if(diffmn.gt.difmax) difmax=diffmn adiff(iuse)=diffmn itag(iuse)=n avedif=avedif+diffmn endif enddo avedif=avedif/iuse ! avedif=(avedif+difmax)/2 emin=9.d190 do i=1,iuse ! print *,"i, adiff(i),avedif : ",i,adiff(i),avedif if(adiff(i).ge.avedif) then itagi=itag(i) benei=bene(itagi) ! print *,"i, benei,emin : ",i,benei,emin if(benei.lt.emin) then emin=benei is(imade1)=itagi endif endif enddo if(ik.eq.0) iuse=iuse-1 ! print *, "exiting select_iseed_min",is(imade1) return end subroutine select_iseed_min !----------------------------------------------------------------------------- subroutine select_iseed_far(imade1,ik) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer :: imade1,ik,n,imade,m real(kind=8) :: dmax,diffmn,diff dmax=-9.d190 do n=1,nbank if(ibank(n).eq.ik) then diffmn=9.d190 do imade=1,imade1-1 ! m=nbank+imade ! call get_diff12(bvar(1,1,1,n),bvar(1,1,1,m),diff,idiff) m=is(imade) diff=dij(n,m) if(diff.lt.diffmn) diffmn=diff enddo endif if(diffmn.gt.dmax) then dmax=diffmn is(imade1)=n endif enddo return end subroutine select_iseed_far !----------------------------------------------------------------------------- subroutine find_min ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer :: i real(kind=8) :: benei ebmin=9.d190 do i=1,nbank benei=bene(i) if(benei.lt.ebmin) then ebmin=benei ibmin=i endif enddo return end subroutine find_min !----------------------------------------------------------------------------- subroutine find_max ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer :: i real(kind=8) :: benei ebmax=-9.d190 do i=1,nbank benei=bene(i) if(benei.gt.ebmax) then ebmax=benei ibmax=i endif enddo return end subroutine find_max !----------------------------------------------------------------------------- subroutine get_diff ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer :: i,i1,i2 real(kind=8) :: tdiff,difmin,diff tdiff=0.d0 difmin=9.d190 do i1=1,nbank-1 do i2=i1+1,nbank if(jbank(i1).eq.0.or.jbank(i2).eq.0) then call get_diff12(bvar(1,1,1,i1),bvar(1,1,1,i2),diff) dij(i1,i2)=diff dij(i2,i1)=diff else diff=dij(i1,i2) endif tdiff=tdiff+diff if(diff.lt.difmin) difmin=diff enddo dij(i1,i1)=0.0 enddo do i=1,nbank jbank(i)=1 enddo avedif=tdiff/nbank/(nbank-1)*2 return end subroutine get_diff !----------------------------------------------------------------------------- subroutine estimate_cutdif(adif,xct,cutdifr) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer :: nexp real(kind=8) :: adif,xct,cutdifr,ctdif1,exponent ctdif1=adif/cut2 exponent = cutdifr*cut1/adif exponent = dlog(exponent)/dlog(xct) nexp=exponent+0.25 cutdif= adif/cut1*xct**nexp if(cutdif.lt.ctdif1) cutdif=ctdif1 return end subroutine estimate_cutdif !----------------------------------------------------------------------------- subroutine get_is_max(idum,ifar,n,imade,k) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' integer :: idum,ifar,n,imade,k,i,j real(kind=8) :: emax do i=imade+1,n emax=-9.d190 do j=1,nbank if(ibank(j).eq.k .and. bene(j).gt.emax) then emax=bene(j) is(i)=j endif enddo call save_is(i) enddo return end subroutine get_is_max !----------------------------------------------------------------------------- ! csa.f !----------------------------------------------------------------------------- subroutine make_array use energy_data, only: itype ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.CSA' integer :: k,j,i,indg !cccccccccccccccccccccccc ! Level-2: group !cccccccccccccccccccccccc indg=0 do k=1,numch !cccccccccccccccccccccccccccccccccccccccc ! Groups the THETAs and the GAMMAs do j=2,nres-1 indg=indg+1 if (j.lt.nres-1) then ngroup(indg)=2 else ngroup(indg)=1 endif do i=1,ngroup(indg) igroup(1,i,indg)=i igroup(2,i,indg)=j igroup(3,i,indg)=k enddo enddo !cccccccccccccccccccccccccccccccccccccccc enddo ! Groups the ALPHAs and the BETAs do k=1,numch do j=2,nres-1 if(itype(j).ne.10) then indg=indg+1 ngroup(indg)=2 do i=1,ngroup(indg) igroup(1,i,indg)=i+2 igroup(2,i,indg)=j igroup(3,i,indg)=k enddo endif enddo enddo ntotgr=indg write(iout,*) write(iout,*) "# of groups: ",ntotgr do i=1,ntotgr write(iout,41) i,ngroup(i),((igroup(k,j,i),k=1,3),j=1,ngroup(i)) enddo ! close(iout) 40 format(i3,3x,3i3) 41 format(2i3,3x,6(3i3,2x)) return end subroutine make_array !----------------------------------------------------------------------------- subroutine make_ranvar(n,m,idum) use geometry_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.BANK' integer :: n,m,j,idum,itrial,jeden ! al m=0 print *,'HOHOHOHO Make_RanVar!!!!!',n,m itrial=0 do while(m.lt.n .and. itrial.le.10000) itrial=itrial+1 jeden=1 call gen_rand_conf(jeden,*10) ! call intout m=m+1 do j=2,nres-1 dihang_in(1,j,1,m)=theta(j+1) dihang_in(2,j,1,m)=phi(j+2) dihang_in(3,j,1,m)=alph(j) dihang_in(4,j,1,m)=omeg(j) enddo dihang_in(2,nres-1,1,m)=0.0d0 goto 20 10 write (iout,*) 'Failed to generate conformation #',m+1,& ' itrial=',itrial 20 continue enddo print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial return end subroutine make_ranvar !----------------------------------------------------------------------------- subroutine make_ranvar_reg(n,idum) use geometry_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.BANK' ! include 'COMMON.GEO' integer :: n,idum,j,m,itrial,jeden m=0 print *,'HOHOHOHO Make_RanVar!!!!!' itrial=0 do while(m.lt.n .and. itrial.le.10000) itrial=itrial+1 jeden=1 call gen_rand_conf(jeden,*10) ! call intout m=m+1 do j=2,nres-1 dihang_in(1,j,1,m)=theta(j+1) dihang_in(2,j,1,m)=phi(j+2) dihang_in(3,j,1,m)=alph(j) dihang_in(4,j,1,m)=omeg(j) if(m.le.n*0.1) then dihang_in(1,j,1,m)=90.0*deg2rad dihang_in(2,j,1,m)=50.0*deg2rad endif enddo dihang_in(2,nres-1,1,m)=0.0d0 goto 20 10 write (iout,*) 'Failed to generate conformation #',m+1,& ' itrial=',itrial 20 continue enddo print *,'Make_RanVar!!!!! m=',m,' itrial=',itrial return end subroutine make_ranvar_reg !----------------------------------------------------------------------------- ! diff12.f !----------------------------------------------------------------------------- subroutine get_diff12(aarray,barray,diff) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' integer :: k,j,i real(kind=8),dimension(mxang,nres,mxch) :: aarray,barray !(mxang,maxres,mxch) real(kind=8) :: diff,dif diff=0.d0 do k=1,numch do j=2,nres-1 ! do i=1,4 ! do i=1,2 do i=1,ndiff dif=rad2deg*dabs(aarray(i,j,k)-barray(i,j,k)) if(dif.gt.180.) dif=360.-dif if (dif.gt.diffcut) diff=diff+dif enddo enddo enddo return end subroutine get_diff12 !----------------------------------------------------------------------------- ! indexx.f !----------------------------------------------------------------------------- subroutine indexx(n,arr,indx) ! implicit real*8 (a-h,o-z) INTEGER :: n,indx(n) REAL(kind=8) :: arr(n) ! PARAMETER (M=7,NSTACK=50) integer,PARAMETER :: M=7,NSTACK=500 INTEGER :: i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK) REAL(kind=8) :: a do 11 j=1,n indx(j)=j 11 continue jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 13 j=l+1,ir indxt=indx(j) a=arr(indxt) do 12 i=j-1,1,-1 if(arr(indx(i)).le.a)goto 2 indx(i+1)=indx(i) 12 continue i=0 2 indx(i+1)=indxt 13 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 itemp=indx(k) indx(k)=indx(l+1) indx(l+1)=itemp if(arr(indx(l+1)).gt.arr(indx(ir)))then itemp=indx(l+1) indx(l+1)=indx(ir) indx(ir)=itemp endif if(arr(indx(l)).gt.arr(indx(ir)))then itemp=indx(l) indx(l)=indx(ir) indx(ir)=itemp endif if(arr(indx(l+1)).gt.arr(indx(l)))then itemp=indx(l+1) indx(l+1)=indx(l) indx(l)=itemp endif i=l+1 j=ir indxt=indx(l) a=arr(indxt) 3 continue i=i+1 if(arr(indx(i)).lt.a)goto 3 4 continue j=j-1 if(arr(indx(j)).gt.a)goto 4 if(j.lt.i)goto 5 itemp=indx(i) indx(i)=indx(j) indx(j)=itemp goto 3 5 indx(l)=indx(j) indx(j)=indxt jstack=jstack+2 if(jstack.gt.NSTACK)pause 'NSTACK too small in indexx' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 end subroutine indexx ! (C) Copr. 1986-92 Numerical Recipes Software *11915aZ%. !----------------------------------------------------------------------------- ! minim_jlee.F !----------------------------------------------------------------------------- subroutine minim_jlee use minim_data use MPI_data use energy_data use compare_data use control_data use geometry_data, only: nvar,nphi use geometry, only:dist use energy, only:fdum use control, only:init_int_table use minimm, only:sumsl,deflt ! controls minimization and sorting routines ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.MINIM' ! include 'COMMON.CONTROL' include 'mpif.h' integer,parameter :: liv=60 integer :: lv ! external func,gradient!,fdum !use minim & energy ! real(kind=4) :: ran1,ran2,ran3 ! include 'COMMON.SETUP' ! include 'COMMON.GEO' ! include 'COMMON.FFIELD' ! include 'COMMON.SBRIDGE' ! include 'COMMON.DISTFIT' ! include 'COMMON.CHAIN' integer,dimension(mpi_status_size) :: muster real(kind=8),dimension(6*nres) :: var !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(mxch*(mxch+1)/2+1) :: erg real(kind=8),dimension(6*nres) :: var2 !(maxvar) (maxvar=6*maxres) integer,dimension(nres) :: iffr !(maxres) integer,dimension((nres-1)*(nres-2)/2) :: ihpbt,jhpbt !(maxdim) (maxdim=(maxres-1)*(maxres-2)/2) real(kind=8),dimension(6*nres) :: d,garbage !(maxvar) (maxvar=6*maxres) !el real(kind=8),dimension(1:lv+1) :: v real(kind=8) :: energia(0:n_ene),time0s,time1s integer,dimension(9) :: indx integer,dimension(12) :: info integer,dimension(liv) :: iv integer :: idum(1) real(kind=8) :: rdum(1) integer,dimension(2,12*nres) :: icont_ !(2,maxcont)(maxcont=12*maxres) logical :: fail !check_var, integer :: iloop(2) !el common /przechowalnia/ v integer :: i,j,ierr,n,nfun,nft_sc,nf,ierror,ierrcode real(kind=8) :: rad,eee,etot !,fdum !el from subroutine parmread ! Define the constants of the disulfide bridge ! Old arbitrary potential real(kind=8),parameter :: dbr=4.20D0 real(kind=8),parameter :: fbr=3.30D0 !----------------- lv=77+(6*nres)*(6*nres+17)/2 !77+maxvar*(maxvar+17)/2 (maxvar=6*maxres) data rad /1.745329252d-2/ ! receive # of start ! print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun, ! & ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf if (.not. allocated(v)) allocate(v(1:lv)) nhpb0=nhpb 10 continue time0s=MPI_WTIME() ! print *, 'MINIM_JLEE: ',me,' is waiting' call mpi_recv(info,12,mpi_integer,king,idint,CG_COMM,& muster,ierr) time1s=MPI_WTIME() write (iout,'(a12,f10.4,a4)')'Waiting for ',time1s-time0s,' sec' call flush(iout) n=info(1) ! print *, 'MINIM_JLEE: ',me,' received: ',n !rc if (ierr.ne.0) go to 100 ! if # = 0, return if (n.eq.0) then write (iout,*) 'Finishing minim_jlee - signal',n,' from master' call flush(iout) return endif nfun=0 IF (n.lt.0) THEN call mpi_recv(var,nvar,mpi_double_precision,& king,idreal,CG_COMM,muster,ierr) call mpi_recv(iffr,nres,mpi_integer,& king,idint,CG_COMM,muster,ierr) call mpi_recv(var2,nvar,mpi_double_precision,& king,idreal,CG_COMM,muster,ierr) ELSE ! receive initial values of variables call mpi_recv(var,nvar,mpi_double_precision,& king,idreal,CG_COMM,muster,ierr) !rc if (ierr.ne.0) go to 100 ENDIF if(vdisulf.and.info(2).ne.-1) then if(info(4).ne.0)then call mpi_recv(ihpbt,info(4),mpi_integer,& king,idint,CG_COMM,muster,ierr) call mpi_recv(jhpbt,info(4),mpi_integer,& king,idint,CG_COMM,muster,ierr) endif endif IF (n.lt.0) THEN n=-n nhpb=nhpb0 link_start=1 link_end=nhpb call init_int_table call contact_cp(var,var2,iffr,nfun,n) ENDIF if(vdisulf.and.info(2).ne.-1) then nss=0 if(info(4).ne.0)then !d write(iout,*) 'SS=',info(4),'N=',info(1),'IT=',info(2) call var_to_geom(nvar,var) call chainbuild do i=1,info(4) if (dist(ihpbt(i),jhpbt(i)).lt.7.0) then nss=nss+1 ihpb(nss)=ihpbt(i) jhpb(nss)=jhpbt(i) !d write(iout,*) 'SS mv=',info(3), !d & ihpb(nss)-nres,jhpb(nss)-nres, !d & dist(ihpb(nss),jhpb(nss)) dhpb(nss)=dbr forcon(nss)=fbr else !d write(iout,*) 'rm SS mv=',info(3), !d & ihpbt(i)-nres,jhpbt(i)-nres,dist(ihpbt(i),jhpbt(i)) endif enddo endif nhpb=nss link_start=1 link_end=nhpb call init_int_table endif if (info(3).eq.14) then write(iout,*) 'calling local_move',info(7),info(8) call local_move_init(.false.) call var_to_geom(nvar,var) call local_move(info(7),info(8),20d0,50d0) call geom_to_var(nvar,var) endif if (info(3).eq.16) then write(iout,*) 'calling beta_slide',info(7),info(8),& info(10), info(11), info(12) call var_to_geom(nvar,var) call beta_slide(info(7),info(8),info(10),info(11),info(12), & nfun,n) call geom_to_var(nvar,var) endif if (info(3).eq.17) then write(iout,*) 'calling beta_zip',info(7),info(8) call var_to_geom(nvar,var) call beta_zip(info(7),info(8),nfun,n) call geom_to_var(nvar,var) endif !rc overlap test if (overlapsc) then call var_to_geom(nvar,var) call chainbuild call etotal(energia) nfun=nfun+1 if (energia(1).eq.1.0d20) then info(3)=-info(3) write (iout,'(a,1pe14.5)')'#OVERLAP evdw=1d20',energia(1) call overlap_sc(fail) if(.not.fail) then call geom_to_var(nvar,var) call etotal(energia) nfun=nfun+1 write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1) else v(10)=1.0d20 iv(1)=-1 goto 201 endif endif endif if (searchsc) then call var_to_geom(nvar,var) call sc_move(2,nres-1,1,10d0,nft_sc,etot) call geom_to_var(nvar,var) !d write(iout,*) 'sc_move',nft_sc,etot endif if (check_var(var,info)) then v(10)=1.0d21 iv(1)=6 goto 201 endif !rc ! write (iout,*) 'MINIM_JLEE: Processor',me,' nvar',nvar ! write (iout,'(8f10.4)') (var(i),i=1,nvar) ! write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar ! write (*,'(8f10.4)') (var(i),i=1,nvar) do i=1,nvar garbage(i)=var(i) enddo call deflt(2,iv,liv,lv,v) ! 12 means fresh start, dont call deflt iv(1)=12 ! max num of fun calls if (maxfun.eq.0) maxfun=500 iv(17)=maxfun ! max num of iterations if (maxmin.eq.0) maxmin=1000 iv(18)=maxmin ! controls output iv(19)=2 ! selects output unit !d iv(21)=iout iv(21)=0 ! 1 means to print out result iv(22)=0 !d iv(22)=1 ! 1 means to print out summary stats iv(23)=0 ! 1 means to print initial x and d iv(24)=0 ! if(me.eq.3.and.n.eq.255) then ! print *,' CHUJ: stoi' ! iv(21)=6 ! iv(22)=1 ! iv(23)=1 ! iv(24)=1 ! endif ! min val for v(radfac) default is 0.1 v(24)=0.1D0 ! max val for v(radfac) default is 4.0 v(25)=2.0D0 ! v(25)=4.0D0 ! check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) ! the sumsl default is 0.1 v(26)=0.1D0 ! false conv if (act fnctn decrease) .lt. v(34) ! the sumsl default is 100*machep v(34)=v(34)/100.0D0 ! absolute convergence if (tolf.eq.0.0D0) tolf=1.0D-4 v(31)=tolf ! relative convergence if (rtolf.eq.0.0D0) rtolf=1.0D-4 v(32)=rtolf ! controls initial step size v(35)=1.0D-1 ! large vals of d correspond to small components of step do i=1,nphi d(i)=1.0D-1 enddo do i=nphi+1,nvar d(i)=1.0D-1 enddo ! minimize energy ! write (iout,*) 'Processor',me,' nvar',nvar ! write (iout,*) 'Variables BEFORE minimization:' ! write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar) ! print *, 'MINIM_JLEE: ',me,' before SUMSL ' call func(nvar,var,nf,eee,idum,rdum,fdum) nfun=nfun+1 if(eee.ge.1.0d20) then ! print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL' ! print *,' energy before SUMSL =',eee ! print *,' aborting local minimization' iv(1)=-1 v(10)=eee go to 201 endif !t time0s=MPI_WTIME() call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum) !t write(iout,*) 'sumsl time=',MPI_WTIME()-time0s,iv(7),v(10) ! print *, 'MINIM_JLEE: ',me,' after SUMSL ' ! find which conformation was returned from sumsl nfun=nfun+iv(7) ! print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf, ! & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32) ! if (iv(1).ne.4 .or. nf.le.1) then ! write (*,*) 'Processor',me,' something bad in SUMSL',iv(1),nf ! write (*,*) 'Initial Variables' ! write (*,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar) ! write (*,*) 'Variables' ! write (*,'(8f10.4)') (rad2deg*var(i),i=1,nvar) ! write (*,*) 'Vector d' ! write (*,'(8f10.4)') (d(i),i=1,nvar) ! write (iout,*) 'Processor',me,' something bad in SUMSL', ! & iv(1),nf ! write (iout,*) 'Initial Variables' ! write (iout,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar) ! write (iout,*) 'Variables' ! write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar) ! write (iout,*) 'Vector d' ! write (iout,'(8f10.4)') (d(i),i=1,nvar) ! endif ! if (nf.lt.iv(6)-1) then ! recalculate intra- and interchain energies ! call func(nvar,var,nf,v(10),iv,v,fdum) ! else if (nf.eq.iv(6)-1) then ! regenerate conformation ! call var_to_geom(nvar,var) ! call chainbuild ! endif ! change origin and axes to standard ECEPP format ! call var_to_geom(nvar,var) ! write (iout,*) 'MINIM_JLEE after minim: Processor',me,' nvar',nvar ! write (iout,'(8f10.4)') (var(i),i=1,nvar) ! write (iout,*) 'Energy:',v(10) ! send back output ! print *, 'MINIM_JLEE: ',me,' minimized: ',n 201 continue indx(1)=n ! return code: 6-gradient 9-number of ftn evaluation, etc indx(2)=iv(1) ! total # of ftn evaluations (for iwf=0, it includes all minimizations). indx(3)=nfun indx(4)=info(2) indx(5)=info(3) indx(6)=nss indx(7)=info(5) indx(8)=info(6) indx(9)=info(9) call mpi_send(indx,9,mpi_integer,king,idint,CG_COMM,& ierr) ! send back energies ! al & cc ! calculate contact order #ifdef CO_BIAS call contact(.false.,ncont,icont_,co) erg(1)=v(10)-1.0d2*co #else erg(1)=v(10) #endif j=1 call mpi_send(erg,j,mpi_double_precision,king,idreal,& CG_COMM,ierr) #ifdef CO_BIAS call mpi_send(co,j,mpi_double_precision,king,idreal,& CG_COMM,ierr) #endif ! send back values of variables call mpi_send(var,nvar,mpi_double_precision,& king,idreal,CG_COMM,ierr) ! print * , 'MINIM_JLEE: Processor',me,' send erg and var ' if(vdisulf.and.info(2).ne.-1.and.nss.ne.0) then !d call intout !d call chainbuild !d call etotal(energia(0)) !d etot=energia(0) !d call enerprint(energia(0)) call mpi_send(ihpb,nss,mpi_integer,& king,idint,CG_COMM,ierr) call mpi_send(jhpb,nss,mpi_integer,& king,idint,CG_COMM,ierr) endif go to 10 100 print *, ' error in receiving message from emperor', me call mpi_abort(mpi_comm_world,ierror,ierrcode) return 200 print *, ' error in sending message to emperor' call mpi_abort(mpi_comm_world,ierror,ierrcode) return 300 print *, ' error in communicating with emperor' call mpi_abort(mpi_comm_world,ierror,ierrcode) return 956 format (' initial energy could not be calculated',41x) 957 format (80x) 965 format (' convergence code ',i2,' # of function calls ',& i4,' # of gradient calls ',i4,10x) 975 format (' energy ',1p,e12.4,' scaled gradient ',e11.3,32x) end subroutine minim_jlee !----------------------------------------------------------------------------- ! newconf.f !----------------------------------------------------------------------------- subroutine make_var(n,idum,iter_csa) use MD_data use energy_data use compare_data use control_data, only: vdisulf use geometry_data use geometry, only: dist include 'mpif.h' ! use random, only: iran_num,ran_number ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.HAIRPIN' ! include 'COMMON.VAR' ! include 'COMMON.DISTFIT' ! include 'COMMON.GEO' ! include 'COMMON.CONTROL' logical :: nicht_getan,nicht_getan1,fail,lfound integer :: nharp,iharp(4,nres/3),nconf_harp integer :: iisucc(mxio) logical :: ifused(mxio) integer :: nhx_seed(nseed),ihx_seed(4,nres/3,nseed) !max_seed integer :: nhx_use(nseed),ihx_use(0:4,nres/3,nseed) integer :: nlx_seed(nseed),ilx_seed(2,nres/3,nseed),& nlx_use(nseed),ilx_use(nres/3,nseed) ! real(kind=4) :: ran1,ran2 integer :: i,j,k,n,idum,iter_csa,iran,index,n7frag,n8frag,n14frag,& n15frag,nbefrag,nlx_tot,iters,i1,i2,i3,ntot_gen,ngen,iih,& ij,jr,iim,nhx_tot,idummy,iter,iif,iig,icheck,ishift,iang,& n8c,ih_start,ih_end,n7c,index2,isize,nsucc,nacc,j1,nran,& ierror,ierrcode real(kind=8) :: d write (iout,*) 'make_var : nseed=',nseed,'ntry=',n index=0 !----------------------------------------- if (n7.gt.0.or.n8.gt.0.or.n9.gt.0.or.n14.gt.0.or.n15.gt.0 & .or.n16.gt.0.or.n17.gt.0.or.n18.gt.0) & call select_frag(n7frag,n8frag,n14frag,& n15frag,nbefrag,iter_csa) !--------------------------------------------------- ! N18 - random perturbation of one phi(=gamma) angle in a loop ! IF (n18.gt.0) THEN nlx_tot=0 do iters=1,nseed i1=is(iters) nlx_seed(iters)=0 do i2=1,n14frag if (lvar_frag(i2,1).eq.i1) then nlx_seed(iters)=nlx_seed(iters)+5 ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2) ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3) ilx_use(nlx_seed(iters),iters)=5 endif enddo nlx_use(iters)=nlx_seed(iters) nlx_tot=nlx_tot+nlx_seed(iters) enddo if (nlx_tot .ge. n18*nseed) then ntot_gen=n18*nseed else ntot_gen=(nlx_tot/nseed)*nseed endif ngen=0 do while (ngen.lt.ntot_gen) do iters=1,nseed iseed=is(iters) if (nlx_use(iters).gt.0) then nicht_getan=.true. do while (nicht_getan) iih=iran_num(1,nlx_seed(iters)) if (ilx_use(iih,iters).gt.0) then nicht_getan=.false. ilx_use(iih,iters)=ilx_use(iih,iters)-1 nlx_use(iters)=nlx_use(iters)-1 endif enddo ngen=ngen+1 index=index+1 movenx(index)=18 parent(1,index)=iseed parent(2,index)=0 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo jr=iran_num(ilx_seed(1,iih,iters),ilx_seed(2,iih,iters)) d=ran_number(-pi,pi) dihang_in(2,jr-2,1,index)=pinorm(dihang_in(2,jr-2,1,index)+d) if (ngen.eq.ntot_gen) goto 145 endif enddo enddo 145 continue ENDIF !----------------------------------------- ! N17 : zip a beta in a seed by forcing one additional p-p contact ! IF (n17.gt.0) THEN nhx_tot=0 do iters=1,nseed i1=is(iters) nhx_seed(iters)=0 nhx_use(iters)=0 do i2=1,nbefrag if (avar_frag(i2,1).eq.i1) then nhx_seed(iters)=nhx_seed(iters)+1 ihx_use(2,nhx_seed(iters),iters)=1 if (avar_frag(i2,5)-avar_frag(i2,3).le.3.and. & avar_frag(i2,2).gt.1.and.avar_frag(i2,4).lt.nres) then ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1 ihx_use(0,nhx_seed(iters),iters)=1 ihx_use(1,nhx_seed(iters),iters)=0 nhx_use(iters)=nhx_use(iters)+1 else if (avar_frag(i2,4).gt.avar_frag(i2,5)) then if (avar_frag(i2,2).gt.1.and. & avar_frag(i2,4).lt.nres) then ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)+1 ihx_use(0,nhx_seed(iters),iters)=1 ihx_use(1,nhx_seed(iters),iters)=0 nhx_use(iters)=nhx_use(iters)+1 endif if (avar_frag(i2,3).lt.nres.and. & avar_frag(i2,5).gt.1) then ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)-1 ihx_use(0,nhx_seed(iters),iters)= & ihx_use(0,nhx_seed(iters),iters)+1 ihx_use(2,nhx_seed(iters),iters)=0 nhx_use(iters)=nhx_use(iters)+1 endif else if (avar_frag(i2,2).gt.1.and. & avar_frag(i2,4).gt.1) then ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,2)-1 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,4)-1 ihx_use(0,nhx_seed(iters),iters)=1 ihx_use(1,nhx_seed(iters),iters)=0 nhx_use(iters)=nhx_use(iters)+1 endif if (avar_frag(i2,3).lt.nres.and. & avar_frag(i2,5).lt.nres) then ihx_seed(1,nhx_seed(iters),iters)=avar_frag(i2,3)+1 ihx_seed(2,nhx_seed(iters),iters)=avar_frag(i2,5)+1 ihx_use(0,nhx_seed(iters),iters)= & ihx_use(0,nhx_seed(iters),iters)+1 ihx_use(2,nhx_seed(iters),iters)=0 nhx_use(iters)=nhx_use(iters)+1 endif endif endif endif enddo nhx_tot=nhx_tot+nhx_use(iters) !d write (iout,*) "debug N17",iters,nhx_seed(iters), !d & nhx_use(iters),nhx_tot enddo if (nhx_tot .ge. n17*nseed) then ntot_gen=n17*nseed else if (nhx_tot .ge. nseed) then ntot_gen=(nhx_tot/nseed)*nseed else ntot_gen=nhx_tot endif !d write (iout,*) "debug N17==",ntot_gen,nhx_tot,nseed ngen=0 do while (ngen.lt.ntot_gen) do iters=1,nseed iseed=is(iters) if (nhx_use(iters).gt.0) then !d write (iout,*) "debug N17",nhx_use(iters),ngen,ntot_gen !d write (iout,*) "debugN17^", !d & (ihx_use(0,k,iters),k=1,nhx_use(iters)) nicht_getan=.true. do while (nicht_getan) iih=iran_num(1,nhx_seed(iters)) !d write (iout,*) "debugN17^",iih if (ihx_use(0,iih,iters).gt.0) then iim=iran_num(1,2) !d write (iout,*) "debugN17=",iih,nhx_seed(iters) !d write (iout,*) "debugN17-",iim,'##', !d & (ihx_use(k,iih,iters),k=0,2) !d call flush(iout) do while (ihx_use(iim,iih,iters).eq.1) iim=iran_num(1,2) !d write (iout,*) "debugN17-",iim,'##', !d & (ihx_use(k,iih,iters),k=0,2) !d call flush(iout) enddo nicht_getan=.false. ihx_use(iim,iih,iters)=1 ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 nhx_use(iters)=nhx_use(iters)-1 endif enddo ngen=ngen+1 index=index+1 movenx(index)=17 parent(1,index)=iseed parent(2,index)=0 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo if (iim.eq.1) then idata(1,index)=ihx_seed(1,iih,iters) idata(2,index)=ihx_seed(2,iih,iters) else idata(1,index)=ihx_seed(3,iih,iters) idata(2,index)=ihx_seed(4,iih,iters) endif if (ngen.eq.ntot_gen) goto 115 endif enddo enddo 115 continue write (iout,*) "N17",n17," ngen/nseed",ngen/nseed,& ngen,nseed ENDIF !----------------------------------------- ! N16 : slide non local beta in a seed by +/- 1 or +/- 2 ! IF (n16.gt.0) THEN nhx_tot=0 do iters=1,nseed i1=is(iters) nhx_seed(iters)=0 do i2=1,n7frag if (bvar_frag(i2,1).eq.i1) then nhx_seed(iters)=nhx_seed(iters)+1 ihx_seed(1,nhx_seed(iters),iters)=bvar_frag(i2,3) ihx_seed(2,nhx_seed(iters),iters)=bvar_frag(i2,4) ihx_seed(3,nhx_seed(iters),iters)=bvar_frag(i2,5) ihx_seed(4,nhx_seed(iters),iters)=bvar_frag(i2,6) ihx_use(0,nhx_seed(iters),iters)=4 do i3=1,4 ihx_use(i3,nhx_seed(iters),iters)=0 enddo endif enddo nhx_use(iters)=4*nhx_seed(iters) nhx_tot=nhx_tot+nhx_seed(iters) !d write (iout,*) "debug N16",iters,nhx_seed(iters) enddo if (4*nhx_tot .ge. n16*nseed) then ntot_gen=n16*nseed else if (4*nhx_tot .ge. nseed) then ntot_gen=(4*nhx_tot/nseed)*nseed else ntot_gen=4*nhx_tot endif write (iout,*) "debug N16",ntot_gen,4*nhx_tot,nseed ngen=0 do while (ngen.lt.ntot_gen) do iters=1,nseed iseed=is(iters) if (nhx_use(iters).gt.0) then nicht_getan=.true. do while (nicht_getan) iih=iran_num(1,nhx_seed(iters)) if (ihx_use(0,iih,iters).gt.0) then iim=iran_num(1,4) do while (ihx_use(iim,iih,iters).eq.1) !d write (iout,*) iim, !d & ihx_use(0,iih,iters),ihx_use(iim,iih,iters) iim=iran_num(1,4) enddo nicht_getan=.false. ihx_use(iim,iih,iters)=1 ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 nhx_use(iters)=nhx_use(iters)-1 endif enddo ngen=ngen+1 index=index+1 movenx(index)=16 parent(1,index)=iseed parent(2,index)=0 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo do i=1,4 idata(i,index)=ihx_seed(i,iih,iters) enddo idata(5,index)=iim if (ngen.eq.ntot_gen) goto 116 endif enddo enddo 116 continue write (iout,*) "N16",n16," ngen/nseed",ngen/nseed,& ngen,nseed ENDIF !----------------------------------------- ! N15 : copy two 2nd structure elements from 1 or 2 conf. in bank to a seed ! IF (n15.gt.0) THEN do iters=1,nseed iseed=is(iters) do i=1,mxio ifused(i)=.false. enddo do idummy=1,n15 iter=0 84 continue iran=0 iif=iran_num(1,n15frag) do while( (ifused(iif) .or. svar_frag(iif,1).eq.iseed) .and. & iran.le.mxio ) iif=iran_num(1,n15frag) iran=iran+1 enddo if(iran.ge.mxio) goto 811 iran=0 iig=iran_num(1,n15frag) do while( (ifused(iig) .or. svar_frag(iig,1).eq.iseed .or. & .not.(svar_frag(iif,3).lt.svar_frag(iig,2).or. & svar_frag(iig,3).lt.svar_frag(iif,2)) ) .and. & iran.le.mxio ) iig=iran_num(1,n15frag) iran=iran+1 enddo if(iran.ge.mxio) goto 811 index=index+1 movenx(index)=15 parent(1,index)=iseed parent(2,index)=svar_frag(iif,1) parent(3,index)=svar_frag(iig,1) if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif ifused(iif)=.true. ifused(iig)=.true. call newconf_copy(idum,dihang_in(1,1,1,index),& svar_frag(iif,1),svar_frag(iif,2),svar_frag(iif,3)) do j=svar_frag(iig,2),svar_frag(iig,3) do i=1,4 dihang_in(i,j,1,index)=bvar(i,j,1,svar_frag(iig,1)) enddo enddo if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) then index=index-1 ifused(iif)=.false. goto 84 endif endif 811 continue enddo enddo ENDIF !----------------------------------------- ! N14 local_move (Maurizio) for loops in a seed ! IF (n14.gt.0) THEN nlx_tot=0 do iters=1,nseed i1=is(iters) nlx_seed(iters)=0 do i2=1,n14frag if (lvar_frag(i2,1).eq.i1) then nlx_seed(iters)=nlx_seed(iters)+3 ilx_seed(1,nlx_seed(iters),iters)=lvar_frag(i2,2) ilx_seed(2,nlx_seed(iters),iters)=lvar_frag(i2,3) ilx_use(nlx_seed(iters),iters)=3 endif enddo nlx_use(iters)=nlx_seed(iters) nlx_tot=nlx_tot+nlx_seed(iters) !d write (iout,*) "debug N14",iters,nlx_seed(iters) enddo if (nlx_tot .ge. n14*nseed) then ntot_gen=n14*nseed else ntot_gen=(nlx_tot/nseed)*nseed endif !d write (iout,*) "debug N14",ntot_gen,n14frag,nseed ngen=0 do while (ngen.lt.ntot_gen) do iters=1,nseed iseed=is(iters) if (nlx_use(iters).gt.0) then nicht_getan=.true. do while (nicht_getan) iih=iran_num(1,nlx_seed(iters)) if (ilx_use(iih,iters).gt.0) then nicht_getan=.false. ilx_use(iih,iters)=ilx_use(iih,iters)-1 nlx_use(iters)=nlx_use(iters)-1 endif enddo ngen=ngen+1 index=index+1 movenx(index)=14 parent(1,index)=iseed parent(2,index)=0 idata(1,index)=ilx_seed(1,iih,iters) idata(2,index)=ilx_seed(2,iih,iters) if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo if (ngen.eq.ntot_gen) goto 131 endif enddo enddo 131 continue !d write (iout,*) "N14",n14," ngen/nseed",ngen/nseed, !d & ngen,nseed ENDIF !----------------------------------------- ! N9 : shift a helix in a seed ! IF (n9.gt.0) THEN nhx_tot=0 do iters=1,nseed i1=is(iters) nhx_seed(iters)=0 do i2=1,n8frag if (hvar_frag(i2,1).eq.i1) then nhx_seed(iters)=nhx_seed(iters)+1 ihx_seed(1,nhx_seed(iters),iters)=hvar_frag(i2,2) ihx_seed(2,nhx_seed(iters),iters)=hvar_frag(i2,3) ihx_use(0,nhx_seed(iters),iters)=4 do i3=1,4 ihx_use(i3,nhx_seed(iters),iters)=0 enddo endif enddo nhx_use(iters)=4*nhx_seed(iters) nhx_tot=nhx_tot+nhx_seed(iters) !d write (iout,*) "debug N9",iters,nhx_seed(iters) enddo if (4*nhx_tot .ge. n9*nseed) then ntot_gen=n9*nseed else ntot_gen=(4*nhx_tot/nseed)*nseed endif !d write (iout,*) "debug N9",ntot_gen,n8frag,nseed ngen=0 do while (ngen.lt.ntot_gen) do iters=1,nseed iseed=is(iters) if (nhx_use(iters).gt.0) then nicht_getan=.true. do while (nicht_getan) iih=iran_num(1,nhx_seed(iters)) if (ihx_use(0,iih,iters).gt.0) then iim=iran_num(1,4) do while (ihx_use(iim,iih,iters).eq.1) !d write (iout,*) iim, !d & ihx_use(0,iih,iters),ihx_use(iim,iih,iters) iim=iran_num(1,4) enddo nicht_getan=.false. ihx_use(iim,iih,iters)=1 ihx_use(0,iih,iters)=ihx_use(0,iih,iters)-1 nhx_use(iters)=nhx_use(iters)-1 endif enddo ngen=ngen+1 index=index+1 movenx(index)=9 parent(1,index)=iseed parent(2,index)=0 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo jstart=max(nnt,ihx_seed(1,iih,iters)+1) jend=min(nct,ihx_seed(2,iih,iters)) !d write (iout,*) "debug N9",iters,iih,jstart,jend if (iim.eq.1) then ishift=-2 else if (iim.eq.2) then ishift=-1 else if (iim.eq.3) then ishift=1 else if (iim.eq.4) then ishift=2 else write (iout,*) 'CHUJ NASTAPIL: iim=',iim #ifdef MPI !el call mpi_abort(mpi_comm_world,ierror,ierrcode) #endif endif do j=jstart,jend if (itype(j).eq.10) then iang=2 else iang=4 endif do i=1,iang if (j+ishift.ge.nnt.and.j+ishift.le.nct) & dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed) enddo enddo if (ishift.gt.0) then do j=0,ishift-1 if (itype(jend+j).eq.10) then iang=2 else iang=4 endif do i=1,iang if (jend+j.ge.nnt.and.jend+j.le.nct) & dihang_in(i,jstart+j,1,index)=bvar(i,jend+j,1,iseed) enddo enddo else do j=0,-ishift-1 if (itype(jstart+j).eq.10) then iang=2 else iang=4 endif do i=1,iang if (jend+j.ge.nnt.and.jend+j.le.nct) & dihang_in(i,jend+j,1,index)=bvar(i,jstart+j,1,iseed) enddo enddo endif if (ngen.eq.ntot_gen) goto 133 endif enddo enddo 133 continue !d write (iout,*) "N9",n9," ngen/nseed",ngen/nseed, !d & ngen,nseed ENDIF !----------------------------------------- ! N8 : copy a helix from bank to seed ! if (n8.gt.0) then if (n8frag.lt.n8) then write (iout,*) "N8: only ",n8frag,'helices' n8c=n8frag else n8c=n8 endif do iters=1,nseed iseed=is(iters) do i=1,mxio ifused(i)=.false. enddo do idummy=1,n8c iter=0 94 continue iran=0 iif=iran_num(1,n8frag) do while( (ifused(iif) .or. hvar_frag(iif,1).eq.iseed) .and. & iran.le.mxio ) iif=iran_num(1,n8frag) iran=iran+1 enddo if(iran.ge.mxio) goto 911 index=index+1 movenx(index)=8 parent(1,index)=iseed parent(2,index)=hvar_frag(iif,1) if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif ifused(iif)=.true. if (hvar_frag(iif,3)-hvar_frag(iif,2).le.6) then call newconf_copy(idum,dihang_in(1,1,1,index),& hvar_frag(iif,1),hvar_frag(iif,2),hvar_frag(iif,3)) else ih_start=iran_num(hvar_frag(iif,2),hvar_frag(iif,3)-6) ih_end=iran_num(ih_start,hvar_frag(iif,3)) call newconf_copy(idum,dihang_in(1,1,1,index),& hvar_frag(iif,1),ih_start,ih_end) endif iter=iter+1 if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) then index=index-1 ifused(iif)=.false. goto 94 endif endif 911 continue enddo enddo endif !----------------------------------------- ! N7 : copy nonlocal beta fragment from bank to seed ! if (n7.gt.0) then if (n7frag.lt.n7) then write (iout,*) "N7: only ",n7frag,'nonlocal fragments' n7c=n7frag else n7c=n7 endif do i=1,nres do j=1,mxio2 iff_in(i,j)=0 enddo enddo index2=0 do i=1,mxio isend2(i)=0 enddo do iters=1,nseed iseed=is(iters) do i=1,mxio ifused(i)=.false. enddo do idummy=1,n7c iran=0 iif=iran_num(1,n7frag) do while( (ifused(iif) .or. bvar_frag(iif,1).eq.iseed) .and. & iran.le.mxio ) iif=iran_num(1,n7frag) iran=iran+1 enddo !d write (*,'(3i5,l,4i5)'),iters,idummy,iif,ifused(iif), !d & bvar_frag(iif,1),iseed,iran,index2 if(iran.ge.mxio) goto 999 if(index2.ge.mxio2) goto 999 index=index+1 movenx(index)=7 parent(1,index)=iseed parent(2,index)=bvar_frag(iif,1) index2=index2+1 isend2(index)=index2 ifused(iif)=.true. if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in2(i,j,k,index2)=bvar(i,j,k,bvar_frag(iif,1)) enddo enddo enddo if (bvar_frag(iif,2).eq.4) then do i=bvar_frag(iif,3),bvar_frag(iif,4) iff_in(i,index2)=1 enddo if (bvar_frag(iif,5).lt.bvar_frag(iif,6)) then !d print *,'###',bvar_frag(iif,3),bvar_frag(iif,4), !d & bvar_frag(iif,5),bvar_frag(iif,6) do i=bvar_frag(iif,5),bvar_frag(iif,6) iff_in(i,index2)=1 enddo else !d print *,'###',bvar_frag(iif,3),bvar_frag(iif,4), !d & bvar_frag(iif,6),bvar_frag(iif,5) do i=bvar_frag(iif,6),bvar_frag(iif,5) iff_in(i,index2)=1 enddo endif endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo 999 continue enddo enddo endif !----------------------------------------------- ! N6 : copy random continues fragment from bank to seed ! do iters=1,nseed iseed=is(iters) do idummy=1,n6 isize=(is2-is1+1)*ran1(idum)+is1 index=index+1 movenx(index)=6 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif iter=0 104 continue if(icycle.le.0) then i1=nconf* ran1(idum)+1 i1=nbank-nconf+i1 else i1=nbank* ran1(idum)+1 endif if(i1.eq.iseed) goto 104 iter=iter+1 call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize) parent(1,index)=iseed parent(2,index)=i1 if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) goto 104 endif enddo enddo !----------------------------------------- if (n3.gt.0.or.n4.gt.0) call gen_hairpin nconf_harp=0 do iters=1,nseed if (nharp_seed(iters).gt.0) nconf_harp=nconf_harp+1 enddo !----------------------------------------- ! N3 : copy hairpin from bank to seed ! do iters=1,nseed iseed=is(iters) nsucc=0 nacc=0 do idummy=1,n3 index=index+1 iter=0 124 continue if(icycle.le.0) then i1=nconf* ran1(idum)+1 i1=nbank-nconf+i1 else i1=nbank* ran1(idum)+1 endif if(i1.eq.iseed) goto 124 do k=1,nsucc if (i1.eq.iisucc(k).and.nsucc.lt.nconf_harp-1) goto 124 enddo nsucc=nsucc+1 iisucc(nsucc)=i1 iter=iter+1 call newconf_residue_hairpin(idum,dihang_in(1,1,1,index),& i1,fail) if (fail) then if (icycle.le.0 .and. nsucc.eq.nconf .or. & icycle.gt.0 .and. nsucc.eq.nbank) then index=index-1 goto 125 else goto 124 endif endif if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) goto 124 endif movenx(index)=3 parent(1,index)=iseed parent(2,index)=i1 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif nacc=nacc+1 enddo ! if not enough hairpins, supplement with windows 125 continue !dd if (n3.ne.0) write (iout,*) "N3",n3," nsucc",nsucc," nacc",nacc do idummy=nacc+1,n3 isize=(is2-is1+1)*ran1(idum)+is1 index=index+1 movenx(index)=6 parent(1,index)=iseed parent(2,index)=i1 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif iter=0 114 continue if(icycle.le.0) then i1=nconf* ran1(idum)+1 i1=nbank-nconf+i1 else i1=nbank* ran1(idum)+1 endif if(i1.eq.iseed) goto 114 iter=iter+1 call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize) if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) goto 114 endif enddo enddo !----------------------------------------- ! N4 : shift a turn in hairpin in seed ! IF (N4.GT.0) THEN if (4*nharp_tot .ge. n4*nseed) then ntot_gen=n4*nseed else ntot_gen=(4*nharp_tot/nseed)*nseed endif ngen=0 do while (ngen.lt.ntot_gen) do iters=1,nseed iseed=is(iters) ! write (iout,*) 'iters',iters,' iseed',iseed,' nharp_seed', ! & nharp_seed(iters),' nharp_use',nharp_use(iters), ! & ' ntot_gen',ntot_gen ! write (iout,*) 'iharp_use(0)', ! & (iharp_use(0,k,iters),k=1,nharp_seed(iters)) if (nharp_use(iters).gt.0) then nicht_getan=.true. do while (nicht_getan) iih=iran_num(1,nharp_seed(iters)) ! write (iout,*) 'iih',iih,' iharp_use', ! & (iharp_use(k,iih,iters),k=1,4) if (iharp_use(0,iih,iters).gt.0) then nicht_getan1=.true. do while (nicht_getan1) iim=iran_num(1,4) nicht_getan1=iharp_use(iim,iih,iters).eq.1 enddo nicht_getan=.false. iharp_use(iim,iih,iters)=1 iharp_use(0,iih,iters)=iharp_use(0,iih,iters)-1 nharp_use(iters)=nharp_use(iters)-1 !dd write (iout,'(a16,i3,a5,i2,a10,2i4)') !dd & 'N4 selected hairpin',iih,' move',iim,' iharp_seed', !dd & iharp_seed(1,iih,iters),iharp_seed(2,iih,iters) endif enddo ngen=ngen+1 index=index+1 movenx(index)=4 parent(1,index)=iseed parent(2,index)=0 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo jstart=iharp_seed(1,iih,iters)+1 jend=iharp_seed(2,iih,iters) if (iim.eq.1) then ishift=-2 else if (iim.eq.2) then ishift=-1 else if (iim.eq.3) then ishift=1 else if (iim.eq.4) then ishift=2 else write (iout,*) 'CHUJ NASTAPIL: iim=',iim #ifdef MPI !el call mpi_abort(mpi_comm_world,ierror,ierrcode) #endif !el endif ! write (iout,*) 'jstart',jstart,' jend',jend,' ishift',ishift ! write (iout,*) 'Before turn shift' ! do j=2,nres-1 ! theta(j+1)=dihang_in(1,j,1,index) ! phi(j+2)=dihang_in(2,j,1,index) ! alph(j)=dihang_in(3,j,1,index) ! omeg(j)=dihang_in(4,j,1,index) ! enddo ! call intout do j=jstart,jend if (itype(j).eq.10) then iang=2 else iang=4 endif do i=1,iang if (j+ishift.ge.nnt.and.j+ishift.le.nct) & dihang_in(i,j+ishift,1,index)=bvar(i,j,1,iseed) enddo enddo ! write (iout,*) 'After turn shift' ! do j=2,nres-1 ! theta(j+1)=dihang_in(1,j,1,index) ! phi(j+2)=dihang_in(2,j,1,index) ! alph(j)=dihang_in(3,j,1,index) ! omeg(j)=dihang_in(4,j,1,index) ! enddo ! call intout if (ngen.eq.ntot_gen) goto 135 endif enddo enddo ! if not enough hairpins, supplement with windows ! write (iout,*) 'end of enddo' 135 continue !dd write (iout,*) "N4",n4," ngen/nseed",ngen/nseed, !dd & ngen,nseed do iters=1,nseed iseed=is(iters) do idummy=ngen/nseed+1,n4 isize=(is2-is1+1)*ran1(idum)+is1 index=index+1 movenx(index)=6 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif iter=0 134 continue if(icycle.le.0) then i1=nconf* ran1(idum)+1 i1=nbank-nconf+i1 else i1=nbank* ran1(idum)+1 endif if(i1.eq.iseed) goto 134 iter=iter+1 call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize) parent(1,index)=iseed parent(2,index)=i1 if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) goto 134 endif enddo enddo ENDIF !----------------------------------------- ! N5 : copy one residue from bank to seed (normally switched off - use N1) ! do iters=1,nseed iseed=is(iters) isize=1 do i=1,n5 index=index+1 movenx(index)=5 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif iter=0 105 continue if(icycle.le.0) then i1=nconf* ran1(idum)+1 i1=nbank-nconf+i1 else i1=nbank* ran1(idum)+1 endif if(i1.eq.iseed) goto 105 iter=iter+1 call newconf_residue(idum,dihang_in(1,1,1,index),i1,isize) parent(1,index)=iseed parent(2,index)=i1 if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) goto 105 endif enddo enddo !----------------------------------------- ! N2 : copy backbone of one residue from bank or first bank to seed ! (normally switched off - use N1) ! do iters=1,nseed iseed=is(iters) do i=n2,1,-1 if(icycle.le.0.and.iuse.gt.nconf-irr) then iseed=ran1(idum)*nconf+1 iseed=nbank-nconf+iseed endif index=index+1 movenx(index)=2 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif iter=0 102 i1= ran1(idum)*nbank+1 if(i1.eq.iseed) goto 102 iter=iter+1 if(icycle.le.0.and.iuse.gt.nconf-irr) then nran=mod(i-1,nran0)+3 call newconf1arr(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=-iseed parent(2,index)=-i1 else if(icycle.le.0.and.iters.le.iuse) then nran=mod(i-1,nran0)+1 call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=iseed parent(2,index)=-i1 else nran=mod(i-1,nran1)+1 if(ran1(idum).lt.0.5) then call newconf1abr(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=iseed parent(2,index)=-i1 else call newconf1abb(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=iseed parent(2,index)=i1 endif endif if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) goto 102 endif enddo enddo !----------------------------------------- ! N1 : copy backbone or sidechain of one residue from bank or ! first bank to seed ! do iters=1,nseed iseed=is(iters) do i=n1,1,-1 if(icycle.le.0.and.iuse.gt.nconf-irr) then iseed=ran1(idum)*nconf+1 iseed=nbank-nconf+iseed endif index=index+1 movenx(index)=1 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif iter=0 101 i1= ran1(idum)*nbank+1 if(i1.eq.iseed) goto 101 iter=iter+1 if(icycle.le.0.and.iuse.gt.nconf-irr) then nran=mod(i-1,nran0)+3 call newconf1rr(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=-iseed parent(2,index)=-i1 else if(icycle.le.0.and.iters.le.iuse) then nran=mod(i-1,nran0)+1 call newconf1br(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=iseed parent(2,index)=-i1 else nran=mod(i-1,nran1)+1 if(ran1(idum).lt.0.5) then call newconf1br(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=iseed parent(2,index)=-i1 else call newconf1bb(idum,dihang_in(1,1,1,index),nran,i1) parent(1,index)=iseed parent(2,index)=i1 endif endif if(iter.lt.10) then call check_old(icheck,index) if(icheck.eq.1) goto 101 endif enddo enddo !----------------------------------------- ! N0 just all seeds ! IF (n0.gt.0) THEN do iters=1,nseed iseed=is(iters) index=index+1 movenx(index)=0 parent(1,index)=iseed parent(2,index)=0 if (vdisulf) then nss_in(index)=bvar_nss(iseed) do ij=1,nss_in(index) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo endif do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo enddo ENDIF !----------------------------------------- if (vdisulf) then do iters=1,nseed iseed=is(iters) do k=1,numch do j=2,nres-1 theta(j+1)=bvar(1,j,k,iseed) phi(j+2)=bvar(2,j,k,iseed) alph(j)=bvar(3,j,k,iseed) omeg(j)=bvar(4,j,k,iseed) enddo enddo call chainbuild !d write(iout,*) 'makevar DYNSS',iseed,'#',bvar_ns(iseed), !d & (bvar_s(k,iseed),k=1,bvar_ns(iseed)), !d & bvar_nss(iseed), !d & (bvar_ss(1,k,iseed)-nres,'-', !d & bvar_ss(2,k,iseed)-nres,k=1,bvar_nss(iseed)) do i1=1,bvar_ns(iseed) ! ! N10 fussion of free halfcysteines in seed ! first select CYS with distance < 7A ! do j1=i1+1,bvar_ns(iseed) if (dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres) & .lt.7.0.and. & iabs(bvar_s(i1,iseed)-bvar_s(j1,iseed)).gt.3) then index=index+1 movenx(index)=10 parent(1,index)=iseed parent(2,index)=0 do ij=1,bvar_nss(iseed) iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) enddo ij=bvar_nss(iseed)+1 nss_in(index)=ij iss_in(ij,index)=bvar_s(i1,iseed)+nres jss_in(ij,index)=bvar_s(j1,iseed)+nres !d write(iout,*) 'makevar NSS0',index, !d & dist(bvar_s(i1,iseed)+nres,bvar_s(j1,iseed)+nres), !d & nss_in(index),iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo endif enddo ! ! N11 type I transdisulfidation ! do j1=1,bvar_nss(iseed) if (dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)) & .lt.7.0.and. & iabs(bvar_s(i1,iseed)-(bvar_ss(1,j1,iseed)-nres)) & .gt.3) then index=index+1 movenx(index)=11 parent(1,index)=iseed parent(2,index)=0 do ij=1,bvar_nss(iseed) if (ij.ne.j1) then iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) endif enddo nss_in(index)=bvar_nss(iseed) iss_in(j1,index)=bvar_s(i1,iseed)+nres jss_in(j1,index)=bvar_ss(1,j1,iseed) if (iss_in(j1,index).gt.jss_in(j1,index)) then iss_in(j1,index)=bvar_ss(1,j1,iseed) jss_in(j1,index)=bvar_s(i1,iseed)+nres endif !d write(iout,*) 'makevar NSS1 #1',index, !d & bvar_s(i1,iseed),bvar_ss(1,j1,iseed)-nres, !d & dist(bvar_s(i1,iseed)+nres,bvar_ss(1,j1,iseed)), !d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres, !d & ij=1,nss_in(index)) do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo endif if (dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)) & .lt.7.0.and. & iabs(bvar_s(i1,iseed)-(bvar_ss(2,j1,iseed)-nres)) & .gt.3) then index=index+1 movenx(index)=11 parent(1,index)=iseed parent(2,index)=0 do ij=1,bvar_nss(iseed) if (ij.ne.j1) then iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) endif enddo nss_in(index)=bvar_nss(iseed) iss_in(j1,index)=bvar_s(i1,iseed)+nres jss_in(j1,index)=bvar_ss(2,j1,iseed) if (iss_in(j1,index).gt.jss_in(j1,index)) then iss_in(j1,index)=bvar_ss(2,j1,iseed) jss_in(j1,index)=bvar_s(i1,iseed)+nres endif !d write(iout,*) 'makevar NSS1 #2',index, !d & bvar_s(i1,iseed),bvar_ss(2,j1,iseed)-nres, !d & dist(bvar_s(i1,iseed)+nres,bvar_ss(2,j1,iseed)), !d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres, !d & ij=1,nss_in(index)) do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo endif enddo enddo ! ! N12 type II transdisulfidation ! do i1=1,bvar_nss(iseed) do j1=i1+1,bvar_nss(iseed) if (dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)) & .lt.7.0.and. & dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)) & .lt.7.0.and. & iabs(bvar_ss(1,i1,iseed)-bvar_ss(1,j1,iseed)) & .gt.3.and. & iabs(bvar_ss(2,i1,iseed)-bvar_ss(2,j1,iseed)) & .gt.3) then index=index+1 movenx(index)=12 parent(1,index)=iseed parent(2,index)=0 do ij=1,bvar_nss(iseed) if (ij.ne.i1 .and. ij.ne.j1) then iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) endif enddo nss_in(index)=bvar_nss(iseed) iss_in(i1,index)=bvar_ss(1,i1,iseed) jss_in(i1,index)=bvar_ss(1,j1,iseed) if (iss_in(i1,index).gt.jss_in(i1,index)) then iss_in(i1,index)=bvar_ss(1,j1,iseed) jss_in(i1,index)=bvar_ss(1,i1,iseed) endif iss_in(j1,index)=bvar_ss(2,i1,iseed) jss_in(j1,index)=bvar_ss(2,j1,iseed) if (iss_in(j1,index).gt.jss_in(j1,index)) then iss_in(j1,index)=bvar_ss(2,j1,iseed) jss_in(j1,index)=bvar_ss(2,i1,iseed) endif !d write(iout,*) 'makevar NSS2 #1',index, !d & bvar_ss(1,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres, !d & dist(bvar_ss(1,i1,iseed),bvar_ss(1,j1,iseed)), !d & bvar_ss(2,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres, !d & dist(bvar_ss(2,i1,iseed),bvar_ss(2,j1,iseed)), !d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres, !d & ij=1,nss_in(index)) do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo endif if (dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)) & .lt.7.0.and. & dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)) & .lt.7.0.and. & iabs(bvar_ss(1,i1,iseed)-bvar_ss(2,j1,iseed)) & .gt.3.and. & iabs(bvar_ss(2,i1,iseed)-bvar_ss(1,j1,iseed)) & .gt.3) then index=index+1 movenx(index)=12 parent(1,index)=iseed parent(2,index)=0 do ij=1,bvar_nss(iseed) if (ij.ne.i1 .and. ij.ne.j1) then iss_in(ij,index)=bvar_ss(1,ij,iseed) jss_in(ij,index)=bvar_ss(2,ij,iseed) endif enddo nss_in(index)=bvar_nss(iseed) iss_in(i1,index)=bvar_ss(1,i1,iseed) jss_in(i1,index)=bvar_ss(2,j1,iseed) if (iss_in(i1,index).gt.jss_in(i1,index)) then iss_in(i1,index)=bvar_ss(2,j1,iseed) jss_in(i1,index)=bvar_ss(1,i1,iseed) endif iss_in(j1,index)=bvar_ss(2,i1,iseed) jss_in(j1,index)=bvar_ss(1,j1,iseed) if (iss_in(j1,index).gt.jss_in(j1,index)) then iss_in(j1,index)=bvar_ss(1,j1,iseed) jss_in(j1,index)=bvar_ss(2,i1,iseed) endif !d write(iout,*) 'makevar NSS2 #2',index, !d & bvar_ss(1,i1,iseed)-nres,bvar_ss(2,j1,iseed)-nres, !d & dist(bvar_ss(1,i1,iseed),bvar_ss(2,j1,iseed)), !d & bvar_ss(2,i1,iseed)-nres,bvar_ss(1,j1,iseed)-nres, !d & dist(bvar_ss(2,i1,iseed),bvar_ss(1,j1,iseed)), !d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres, !d & ij=1,nss_in(index)) do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo endif enddo enddo ! ! N13 removal of disulfide bond ! if (bvar_nss(iseed).gt.0) then i1=bvar_nss(iseed)*ran1(idum)+1 index=index+1 movenx(index)=13 parent(1,index)=iseed parent(2,index)=0 ij=0 do j1=1,bvar_nss(iseed) if (j1.ne.i1) then ij=ij+1 iss_in(ij,index)=bvar_ss(1,j1,iseed) jss_in(ij,index)=bvar_ss(2,j1,iseed) endif enddo nss_in(index)=bvar_nss(iseed)-1 !d write(iout,*) 'NSS3',index,i1, !d & bvar_ss(1,i1,iseed)-nres,'=',bvar_ss(2,i1,iseed)-nres,'#', !d & (iss_in(ij,index)-nres,'-',jss_in(ij,index)-nres, !d & ij=1,nss_in(index)) do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,iseed) enddo enddo enddo endif enddo endif !----------------------------------------- if(index.ne.n) write(iout,*)'make_var : ntry=',index n=index !d do ii=1,n !d write (istat,*) "======== ii=",ii," the dihang array" !d do i=1,nres !d write (istat,'(i5,4f15.5)') i,(dihang_in(k,i,1,ii)*rad2deg,k=1,4) !d enddo !d enddo return end subroutine make_var !----------------------------------------------------------------------------- subroutine check_old(icheck,n) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' integer :: icheck,n,i1,i2,m,j,i real(kind=8) :: ctdif,ctdiff,diff,dif data ctdif /10./ data ctdiff /60./ i1=n do i2=1,n-1 diff=0.d0 do m=1,numch do j=2,nres-1 do i=1,4 dif=rad2deg*dabs(dihang_in(i,j,m,i1)-dihang_in(i,j,m,i2)) if(dif.gt.180.0) dif=360.0-dif if(dif.gt.ctdif) goto 100 diff=diff+dif if(diff.gt.ctdiff) goto 100 enddo enddo enddo icheck=1 return 100 continue enddo icheck=0 return end subroutine check_old !----------------------------------------------------------------------------- subroutine newconf1rr(idum,vvar,nran,i1) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind real(kind=8) :: ctdif,dif ctdif=10. do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=rvar(i,j,k,iseed) enddo enddo enddo do index=1,nran iold(index) = 0 enddo number=ntotgr iter=0 do index=1,nran 10 iran= ran1(idum)*number+1 if(iter.gt.number) return iter=iter+1 if(iter.eq.1) goto 11 do ind=1,index-1 if(iran.eq.iold(ind)) goto 10 enddo 11 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo if(iter.gt.number) goto 20 goto 10 20 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) vvar(i,j,k)=rvar(i,j,k,i1) enddo iold(index)=iran enddo return end subroutine newconf1rr !----------------------------------------------------------------------------- subroutine newconf1br(idum,vvar,nran,i1) use energy_data, only: ndih_nconstr,idih_nconstr use control_data, only: i2ndstr ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.TORCNSTR' ! include 'COMMON.CONTROL' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,nran,i1,iran,index,number,iter,juhc,ind real(kind=8) :: ctdif,dif,rtmp ctdif=10. do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=bvar(i,j,k,iseed) enddo enddo enddo do index=1,nran iold(index) = 0 enddo number=ntotgr iter=0 do index=1,nran 10 iran= ran1(idum)*number+1 if(i2ndstr.gt.0) then rtmp=ran1(idum) if(rtmp.le.rdih_bias) then i=0 do j=1,ndih_nconstr if(igroup(2,1,iran).eq.idih_nconstr(j))i=j enddo if(i.eq.0) then juhc=0 4321 juhc=juhc+1 iran= ran1(idum)*number+1 i=0 do j=1,ndih_nconstr if(igroup(2,1,iran).eq.idih_nconstr(j))i=j enddo if(i.eq.0.or.juhc.lt.1000)goto 4321 if(juhc.eq.1000) then print *, 'move 6 : failed to find unconstrained group' write(iout,*) 'move 6 : failed to find unconstrained group' endif endif endif endif if(iter.gt.number) return iter=iter+1 if(iter.eq.1) goto 11 do ind=1,index-1 if(iran.eq.iold(ind)) goto 10 enddo 11 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo if(iter.gt.number) goto 20 goto 10 20 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) vvar(i,j,k)=rvar(i,j,k,i1) enddo iold(index)=iran enddo return end subroutine newconf1br !----------------------------------------------------------------------------- subroutine newconf1bb(idum,vvar,nran,i1) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind real(kind=8) :: ctdif,dif ctdif=10. do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=bvar(i,j,k,iseed) enddo enddo enddo do index=1,nran iold(index) = 0 enddo number=ntotgr iter=0 do index=1,nran 10 iran= ran1(idum)*number+1 if(iter.gt.number) return iter=iter+1 if(iter.eq.1) goto 11 do ind=1,index-1 if(iran.eq.iold(ind)) goto 10 enddo 11 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo if(iter.gt.number) goto 20 goto 10 20 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) vvar(i,j,k)=bvar(i,j,k,i1) enddo iold(index)=iran enddo return end subroutine newconf1bb !----------------------------------------------------------------------------- subroutine newconf1arr(idum,vvar,nran,i1) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind real(kind=8) :: ctdif,dif ctdif=10. do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=rvar(i,j,k,iseed) enddo enddo enddo do index=1,nran iold(index) = 0 enddo number=nres-2 iter=0 do index=1,nran 10 iran= ran1(idum)*number+1 if(iter.gt.number) return iter=iter+1 if(iter.eq.1) goto 11 do ind=1,index-1 if(iran.eq.iold(ind)) goto 10 enddo 11 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo if(iter.gt.number) goto 20 goto 10 20 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) vvar(i,j,k)=rvar(i,j,k,i1) enddo iold(index)=iran enddo return end subroutine newconf1arr !----------------------------------------------------------------------------- subroutine newconf1abr(idum,vvar,nran,i1) use energy_data, only: ndih_nconstr,idih_nconstr use control_data, only: i2ndstr ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.TORCNSTR' ! include 'COMMON.CONTROL' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind real(kind=8) :: ctdif,dif,rtmp ctdif=10. do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=bvar(i,j,k,iseed) enddo enddo enddo do index=1,nran iold(index) = 0 enddo number=nres-2 iter=0 do index=1,nran 10 iran= ran1(idum)*number+1 if(i2ndstr.gt.0) then rtmp=ran1(idum) if(rtmp.le.rdih_bias) then iran=ran1(idum)*ndih_nconstr+1 iran=idih_nconstr(iran) endif endif if(iter.gt.number) return iter=iter+1 if(iter.eq.1) goto 11 do ind=1,index-1 if(iran.eq.iold(ind)) goto 10 enddo 11 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) dif=rad2deg*dabs(vvar(i,j,k)-rvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo if(iter.gt.number) goto 20 goto 10 20 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) vvar(i,j,k)=rvar(i,j,k,i1) enddo iold(index)=iran enddo return end subroutine newconf1abr !----------------------------------------------------------------------------- subroutine newconf1abb(idum,vvar,nran,i1) use energy_data, only: ndih_nconstr,idih_nconstr use control_data, only: i2ndstr ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.TORCNSTR' ! include 'COMMON.CONTROL' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,nran,i1,iran,index,number,iter,ind real(kind=8) :: ctdif,dif,rtmp ctdif=10. do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=bvar(i,j,k,iseed) enddo enddo enddo do index=1,nran iold(index) = 0 enddo number=nres-2 iter=0 do index=1,nran 10 iran= ran1(idum)*number+1 if(i2ndstr.gt.0) then rtmp=ran1(idum) if(rtmp.le.rdih_bias) then iran=ran1(idum)*ndih_nconstr+1 iran=idih_nconstr(iran) endif endif if(iter.gt.number) return iter=iter+1 if(iter.eq.1) goto 11 do ind=1,index-1 if(iran.eq.iold(ind)) goto 10 enddo 11 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo if(iter.gt.number) goto 20 goto 10 20 continue do ind=1,ngroup(iran) i=igroup(1,ind,iran) j=igroup(2,ind,iran) k=igroup(3,ind,iran) vvar(i,j,k)=bvar(i,j,k,i1) enddo iold(index)=iran enddo return end subroutine newconf1abb !----------------------------------------------------------------------------- subroutine newconf_residue(idum,vvar,i1,isize) use energy_data, only: ndih_nconstr,idih_nconstr use control_data, only: i2ndstr use MPI_data include 'mpif.h' ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.TORCNSTR' ! include 'COMMON.CONTROL' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,i1,isize,iran,number,iter,ind,iend,istart,& ierror,ierrcode real(kind=8) :: ctdif,dif,rtmp ctdif=10. if (iseed.gt.mxio .or. iseed.lt.1) then write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed call mpi_abort(mpi_comm_world,ierror,ierrcode) endif do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=bvar(i,j,k,iseed) enddo enddo enddo k=1 number=nres+isize-2 iter=1 10 iran= ran1(idum)*number+1 if(i2ndstr.gt.0) then rtmp=ran1(idum) if(rtmp.le.rdih_bias) then iran=ran1(idum)*ndih_nconstr+1 iran=idih_nconstr(iran) endif endif istart=iran-isize+1 iend=iran if(istart.lt.2) istart=2 if(iend.gt.nres-1) iend=nres-1 if(iter.eq.1) goto 11 do ind=1,iter-1 if(iran.eq.iold(ind)) goto 10 enddo 11 continue do j=istart,iend do i=1,4 dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo enddo iold(iter)=iran iter=iter+1 if(iter.gt.number) goto 20 goto 10 20 continue do j=istart,iend do i=1,4 vvar(i,j,k)=bvar(i,j,k,i1) enddo enddo return end subroutine newconf_residue !----------------------------------------------------------------------------- subroutine newconf_copy(idum,vvar,i1,istart,iend) use MPI_data include 'mpif.h' ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.TORCNSTR' ! include 'COMMON.CONTROL' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: i,j,k,idum,i1,istart,iend,ierror,ierrcode real(kind=8) :: ctdif,dif ctdif=10. if (iseed.gt.mxio .or. iseed.lt.1) then write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed call mpi_abort(mpi_comm_world,ierror,ierrcode) endif do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=bvar(i,j,k,iseed) enddo enddo enddo do j=istart,iend do i=1,4 vvar(i,j,1)=bvar(i,j,1,i1) enddo enddo return end subroutine newconf_copy !----------------------------------------------------------------------------- subroutine newconf_residue_hairpin(idum,vvar,i1,fail) use geometry_data ! use random, only: iran_num use MPI_data use compare, only:hairpin include 'mpif.h' ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! real(kind=4) :: ran1,ran2 real(kind=8),dimension(mxang,nres,mxch) :: vvar !(mxang,maxres,mxch) integer,dimension(ntotal) :: iold integer :: nharp,iharp(4,nres/3),icipa(nres/3) logical :: fail,not_done integer :: idum,i,j,k,i1,iend,istart,iii,n_used,icount,iih,& ierror,ierrcode real(kind=8) :: ctdif,dif ctdif=10. fail=.false. if (iseed.gt.mxio .or. iseed.lt.1) then write (iout,*) 'Dimension ERROR in NEWCONF: ISEED',iseed call mpi_abort(mpi_comm_world,ierror,ierrcode) endif do k=1,numch do j=2,nres-1 do i=1,4 vvar(i,j,k)=bvar(i,j,k,iseed) enddo enddo enddo do k=1,numch do j=2,nres-1 theta(j+1)=bvar(1,j,k,i1) phi(j+2)=bvar(2,j,k,i1) alph(j)=bvar(3,j,k,i1) omeg(j)=bvar(4,j,k,i1) enddo enddo ! call intout call chainbuild call hairpin(.false.,nharp,iharp) if (nharp.eq.0) then fail=.true. return endif n_used=0 DO III=1,NHARP not_done = .true. icount=0 do while (not_done) icount=icount+1 iih=iran_num(1,nharp) do k=1,n_used if (iih.eq.icipa(k)) then iih=0 goto 22 endif enddo not_done=.false. n_used=n_used+1 icipa(n_used)=iih 22 continue not_done = not_done .and. icount.le.nharp enddo if (iih.eq.0) then write (iout,*) "CHUJ NASTAPIL W NEWCONF_RESIDUE_HAIRPIN!!!!" fail=.true. return endif istart=iharp(1,iih)+1 iend=iharp(2,iih) !dd write (iout,*) "newconf_residue_hairpin: iih",iih, !dd & " istart",istart," iend",iend do k=1,numch do j=istart,iend do i=1,4 dif=rad2deg*dabs(vvar(i,j,k)-bvar(i,j,k,i1)) if(dif.gt.180.) dif=360.-dif if(dif.gt.ctdif) goto 20 enddo enddo enddo goto 10 20 continue do k=1,numch do j=istart,iend do i=1,4 vvar(i,j,k)=bvar(i,j,k,i1) enddo enddo enddo ! do j=1,numch ! do l=2,nres-1 ! write (iout,'(4f8.3)') (rad2deg*vvar(i,l,j),i=1,4) ! enddo ! enddo return 10 continue ENDDO fail=.true. return end subroutine newconf_residue_hairpin !----------------------------------------------------------------------------- subroutine gen_hairpin use geometry_data use MD_data use compare, only:hairpin ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.HAIRPIN' integer :: i1,j,k,iters ! write (iout,*) 'Entering GEN_HAIRPIN' do iters=1,nseed i1=is(iters) do k=1,numch do j=2,nres-1 theta(j+1)=bvar(1,j,k,i1) phi(j+2)=bvar(2,j,k,i1) alph(j)=bvar(3,j,k,i1) omeg(j)=bvar(4,j,k,i1) enddo enddo call chainbuild call hairpin(.false.,nharp_seed(iters),iharp_seed(1,1,iters)) enddo nharp_tot=0 do iters=1,nseed nharp_tot=nharp_tot+nharp_seed(iters) nharp_use(iters)=4*nharp_seed(iters) do j=1,nharp_seed(iters) iharp_use(0,j,iters)=4 do k=1,4 iharp_use(k,j,iters)=0 enddo enddo enddo write (iout,*) 'GEN_HAIRPIN: nharp_tot',nharp_tot !dd do i=1,nseed !dd write (iout,*) 'seed',i !dd write (iout,*) 'nharp_seed',nharp_seed(i), !dd & ' nharp_use',nharp_use(i) !d write (iout,*) 'iharp_seed, iharp_use' !d do j=1,nharp_seed(i) !d write (iout,'(7i3)') iharp_seed(1,j,i),iharp_seed(2,j,i), !d & (iharp_use(k,j,i),k=0,4) !d enddo !dd enddo return end subroutine gen_hairpin !----------------------------------------------------------------------------- subroutine select_frag(nn,nh,nl,ns,nb,i_csa) use geometry_data use MD_data use compare_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.HAIRPIN' ! include 'COMMON.DISTFIT' character(len=50) :: linia integer :: isec(nres) integer :: i,j,i1,k,nn,nh,nl,ns,nb,i_csa,nl1,ns1 nn=0 nh=0 nl=0 ns=0 nb=0 !d write (iout,*) 'Entering select_frag' do i1=1,nbank do i=1,nres isec(i)=0 enddo do k=1,numch do j=2,nres-1 theta(j+1)=bvar(1,j,k,i1) phi(j+2)=bvar(2,j,k,i1) alph(j)=bvar(3,j,k,i1) omeg(j)=bvar(4,j,k,i1) enddo enddo call chainbuild !d write (iout,*) ' -- ',i1,' -- ' call secondary2(.false.) ! ! bvar_frag nn==pair of nonlocal strands in beta sheet (loop>4) ! strands > 4 residues; used by N7 and N16 ! do j=1,nbfrag ! !test 09/12/02 bfrag(2,j)-bfrag(1,j).gt.3 ! do i=bfrag(1,j),bfrag(2,j) isec(i)=1 enddo do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j)) isec(i)=1 enddo if ( (bfrag(3,j).lt.bfrag(4,j) .or. & bfrag(4,j)-bfrag(2,j).gt.4) .and. & bfrag(2,j)-bfrag(1,j).gt.4 ) then nn=nn+1 if (bfrag(3,j).lt.bfrag(4,j)) then write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') & "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,& ",",bfrag(3,j)-1,"-",bfrag(4,j)-1 else write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)') & "select",bfrag(1,j)-1,"-",bfrag(2,j)-1,& ",",bfrag(4,j)-1,"-",bfrag(3,j)-1 endif !d call write_pdb(i_csa*1000+nn+nh,linia,0d0) bvar_frag(nn,1)=i1 bvar_frag(nn,2)=4 do i=1,4 bvar_frag(nn,i+2)=bfrag(i,j) enddo endif enddo ! ! hvar_frag nh==helices; used by N8 and N9 ! do j=1,nhfrag do i=hfrag(1,j),hfrag(2,j) isec(i)=2 enddo if ( hfrag(2,j)-hfrag(1,j).gt.4 ) then nh=nh+1 !d write(linia,'(a6,i3,a1,i3)') !d & "select",hfrag(1,j)-1,"-",hfrag(2,j)-1 !d call write_pdb(i_csa*1000+nn+nh,linia,0d0) hvar_frag(nh,1)=i1 hvar_frag(nh,2)=hfrag(1,j) hvar_frag(nh,3)=hfrag(2,j) endif enddo !v write(iout,'(i4,1pe12.4,1x,1000i1)') !v & i1,bene(i1),(isec(i),i=1,nres) !v write(linia,'(i4,1x,1000i1)') !v & i1,(isec(i),i=1,nres) !v call write_pdb(i_csa*1000+i1,linia,bene(i1)) ! ! lvar_frag nl==loops; used by N14 ! i=1 nl1=nl do while (i.lt.nres) if (isec(i).eq.0) then nl=nl+1 lvar_frag(nl,1)=i1 lvar_frag(nl,2)=i i=i+1 do while (isec(i).eq.0.and.i.le.nres) i=i+1 enddo lvar_frag(nl,3)=i-1 if (lvar_frag(nl,3)-lvar_frag(nl,2).lt.1) nl=nl-1 endif i=i+1 enddo !d write(iout,'(4i5)') (i,(lvar_frag(i,ii),ii=1,3),i=nl1+1,nl) ! ! svar_frag ns==an secondary structure element; used by N15 ! i=1 ns1=ns do while (i.lt.nres) if (isec(i).gt.0) then ns=ns+1 svar_frag(ns,1)=i1 svar_frag(ns,2)=i i=i+1 do while (isec(i).gt.0.and.isec(i-1).eq.isec(i) & .and.i.le.nres) i=i+1 enddo svar_frag(ns,3)=i-1 if (svar_frag(ns,3)-svar_frag(ns,2).lt.1) ns=ns-1 endif if (isec(i).eq.0) i=i+1 enddo !d write(iout,'(4i5)') (i,(svar_frag(i,ii),ii=1,3),i=ns1+1,ns) ! ! avar_frag nb==any pair of beta strands; used by N17 ! do j=1,nbfrag nb=nb+1 avar_frag(nb,1)=i1 do i=1,4 avar_frag(nb,i+1)=bfrag(i,j) enddo enddo enddo return end subroutine select_frag !----------------------------------------------------------------------------- ! together.F !----------------------------------------------------------------------------- subroutine together ! feeds tasks for parallel processing use MPI_data use geometry_data use control_data, only: vdisulf use energy_data use io, only:from_int,write_csa_pdb ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! real(kind=4) :: ran1,ran2 ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.TIME1' ! include 'COMMON.SETUP' ! include 'COMMON.VAR' ! include 'COMMON.GEO' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' real(kind=4) :: tcpu real(kind=8) :: time_start,time_start_c,time0f,time0i logical :: ovrtim,sync_iter,timeout,flag,timeout1 integer,dimension(mpi_status_size) :: muster real(kind=8),dimension(0:100) :: t100 integer,dimension(mxio) :: indx real(kind=8),dimension(6*nres) :: xout !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout integer,dimension(9) :: ind real(kind=8),dimension(2) :: cout real(kind=8),parameter :: rad=1.745329252d-2 integer :: i,m,j,jlee,nft,idum,nrmsdb,nrmsdb1,ierr,ierror,ierrcode,& ntrial,ntry,idum2,imax,idumm,nconfr,iconf,mm,k,im,nst,ifar,& iter,irecv,isent,iw_pdb,nft0i,nft00_c,nft00,ifrom,ij,& ireq,ireq2 real(kind=8) :: adif,p_cut,cutdifr,rmsdbc1c,time1i,ctdif1,xctdif,& time2i,tstart,tend1 !ccccccccccccccccccccccccccccccccccccccccccccccc sync_iter=.true. !el nft=0 !el time_start=0.0d0 IF (ME.EQ.KING) THEN time0f=MPI_WTIME() ilastnstep=1 sync_iter=.false. numch=1 nrmsdb=0 nrmsdb1=0 rmsdbc1c=rmsdbc1 nstep=0 call csa_read call make_array if(iref.ne.0) call from_int(1,0,idum) ! To minimize input conformation (bank conformation) ! Output to $mol.reminimized if (irestart.lt.0) then call read_bank(0,nft,cutdifr) if (irestart.lt.-10) then p_cut=nres*4.d0 call prune_bank(p_cut) return endif call reminimize(jlee) return endif if (irestart.eq.0) then call initial_write nbank=nconf ntbank=nconf if (ntbankm.eq.0) ntbank=0 nstep=0 nft=0 do i=1,mxio ibank(i)=0 jbank(i)=0 enddo else call restart_write !!bankt call read_bankt(jlee,nft,cutdifr) call read_bank(jlee,nft,cutdifr) call read_rbank(jlee,adif) if(iref.ne.0) call from_int(1,0,idum) endif nstmax=nstmax+nstep ntrial=n1+n2+n3+n4+n5+n6+n7+n8 ntry=ntrial+1 ntry=ntry*nseed ! ntrial : number of trial conformations per seed. ! ntry : total number of trial conformations including seed conformations. idum2=-123 ! imax=2**31-1 imax=huge(0) ENDIF call mpi_bcast(jend,1,mpi_integer,0,CG_COMM,ierr) !ccccccccccccccccccccccccccccccccccccccc do 300 jlee=1,jend !ccccccccccccccccccccccccccccccccccccccc 331 continue IF (ME.EQ.KING) THEN if(sync_iter) goto 333 idum=- ran2(idum2)*imax if(jlee.lt.jstart) goto 300 ! Restart the random number generator for conformation generation if(irestart.gt.0) then idum2=idum2+nstep if(idum2.le.0) idum2=-idum2+1 idum=- ran2(idum2)*imax endif idumm=idum call vrndst(idumm) open(icsa_seed,file=csa_seed,status="old") write(icsa_seed,*) "jlee : ",jlee close(icsa_seed) call history_append write(icsa_history,*) "number of procs is ",nodes write(icsa_history,*) jlee,idum,idum2 close(icsa_history) !ccccccccccccccccccccccccccccccccccccccccccccccc 333 icycle=0 call history_append write(icsa_history,*) "nbank is ",nbank close(icsa_history) if(irestart.eq.1) goto 111 if(irestart.eq.2) then icycle=0 do i=1,nbank ibank(i)=1 enddo do i=nbank+1,nbank+nconf ibank(i)=0 enddo endif ! start energy minimization nconfr=max0(nconf+nadd,nodes-1) if (sync_iter) nconf_in=0 ! king-emperor - feed input and sort output write (iout,*) "NCONF_IN",nconf_in m=0 if (nconf_in.gt.0) then ! al 7/2/00 - added possibility to read in some of the initial conformations do m=1,nconf_in read (intin,'(i5)',end=11,err=12) iconf 12 continue write (iout,*) "write READ_ANGLES",iconf,m call read_angles(intin,*11) if (iref.eq.0) then mm=m else mm=m+1 endif do j=2,nres-1 dihang_in(1,j,1,mm)=theta(j+1) dihang_in(2,j,1,mm)=phi(j+2) dihang_in(3,j,1,mm)=alph(j) dihang_in(4,j,1,mm)=omeg(j) enddo enddo ! m goto 13 11 write (iout,*) nconf_in," conformations requested, but only",& m-1," found in the angle file." nconf_in=m-1 13 continue m=nconf_in write (iout,*) nconf_in,& " initial conformations have been read in." endif if (iref.eq.0) then if (nconfr.gt.nconf_in) then call make_ranvar(nconfr,m,idum) write (iout,*) nconfr-nconf_in,& " conformations have been generated randomly." endif else nconfr=nconfr*2 call from_int(nconfr,m,idum) ! call from_pdb(nconfr,idum) endif write (iout,*) 'Exitted from make_ranvar nconfr=',nconfr write (*,*) 'Exitted from make_ranvar nconfr=',nconfr do m=1,nconfr write (iout,*) 'Initial conformation',m write(iout,'(8f10.4)') (rad2deg*dihang_in(1,j,1,m),j=2,nres-1) write(iout,'(8f10.4)') (rad2deg*dihang_in(2,j,1,m),j=2,nres-1) write(iout,'(8f10.4)') (rad2deg*dihang_in(3,j,1,m),j=2,nres-1) write(iout,'(8f10.4)') (rad2deg*dihang_in(4,j,1,m),j=2,nres-1) enddo write(iout,*)'Calling FEEDIN NCONF',nconfr time1i=MPI_WTIME() call feedin(nconfr,nft) write (iout,*) ' Time for first bank min.',MPI_WTIME()-time1i call history_append write(icsa_history,*) jlee,nft,nbank write(icsa_history,851) (etot(i),i=1,nconfr) write(icsa_history,850) (rmsn(i),i=1,nconfr) write(icsa_history,850) (pncn(i),i=1,nconfr) write(icsa_history,*) close(icsa_history) ELSE ! To minimize input conformation (bank conformation) ! Output to $mol.reminimized if (irestart.lt.0) then call reminimize(jlee) return endif if (irestart.eq.1) goto 111 ! soldier - perform energy minimization 334 call minim_jlee ENDIF !cccccccccccccccccccccccccccccccccc ! need to syncronize all procs call mpi_barrier(CG_COMM,ierr) if (ierr.ne.0) then print *, ' cannot synchronize MPI' stop endif !cccccccccccccccccccccccccccccccccc IF (ME.EQ.KING) THEN ! print *,"ok after minim" nstep=nstep+nconf if(irestart.eq.2) then nbank=nbank+nconf ! ntbank=ntbank+nconf if(ntbank.gt.ntbankm) ntbank=ntbankm endif ! print *,"ok before indexx" if(iref.eq.0) then call indexx(nconfr,etot,indx) else ! cc/al 7/6/00 do k=1,nconfr indx(k)=k enddo call indexx(nconfr-nconf_in,rmsn(nconf_in+1),indx(nconf_in+1)) do k=nconf_in+1,nconfr indx(k)=indx(k)+nconf_in enddo ! cc/al ! call indexx(nconfr,rmsn,indx) endif ! print *,"ok after indexx" do im=1,nconf m=indx(im) if (m.gt.mxio .or. m.lt.1) then write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,' M',m call mpi_abort(mpi_comm_world,ierror,ierrcode) endif jbank(im+nbank-nconf)=0 bene(im+nbank-nconf)=etot(m) rene(im+nbank-nconf)=etot(m) !!bankt btene(im)=etot(m) ! brmsn(im+nbank-nconf)=rmsn(m) bpncn(im+nbank-nconf)=pncn(m) rrmsn(im+nbank-nconf)=rmsn(m) rpncn(im+nbank-nconf)=pncn(m) if (im+nbank-nconf.gt.mxio .or. im+nbank-nconf.lt.1) then write (iout,*) 'Dimension ERROR in TOGEHER: IM',im,& ' NBANK',nbank,' NCONF',nconf,' IM+NBANK-NCONF',im+nbank-nconf call mpi_abort(mpi_comm_world,ierror,ierrcode) endif do k=1,numch do j=2,nres-1 do i=1,4 bvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m) rvar(i,j,k,im+nbank-nconf)=dihang(i,j,k,m) !!bankt btvar(i,j,k,im)=dihang(i,j,k,m) ! enddo enddo enddo if(iref.eq.1) then if(brmsn(im+nbank-nconf).gt.rmscut.or. & bpncn(im+nbank-nconf).lt.pnccut) ibank(im+nbank-nconf)=9 endif if(vdisulf) then bvar_ns(im+nbank-nconf)=ns-2*nss k=0 do i=1,ns j=1 do while( iss(i).ne.ihpb(j)-nres .and. & iss(i).ne.jhpb(j)-nres .and. j.le.nss) j=j+1 enddo if (j.gt.nss) then k=k+1 bvar_s(k,im+nbank-nconf)=iss(i) endif enddo endif bvar_nss(im+nbank-nconf)=nss do i=1,nss bvar_ss(1,i,im+nbank-nconf)=ihpb(i) bvar_ss(2,i,im+nbank-nconf)=jhpb(i) enddo enddo ENDIF 111 continue IF (ME.EQ.KING) THEN call find_max call find_min call get_diff if(nbank.eq.nconf.and.irestart.eq.0) then adif=avedif endif cutdif=adif/cut1 ctdif1=adif/cut2 !d print *,"adif,xctdif,cutdifr" !d print *,adif,xctdif,cutdifr nst=ntotal/ntrial/nseed xctdif=(cutdif/ctdif1)**(-1.0/nst) if(irestart.ge.1) call estimate_cutdif(adif,xctdif,cutdifr) ! print *,"ok after estimate" irestart=0 call write_rbank(jlee,adif,nft) call write_bank(jlee,nft) !!bankt call write_bankt(jlee,nft) ! call write_bank1(jlee) call history_append write(icsa_history,*) "xctdif: ", xctdif,nst,adif/cut1,ctdif1 write(icsa_history,851) (bene(i),i=1,nbank) write(icsa_history,850) (brmsn(i),i=1,nbank) write(icsa_history,850) (bpncn(i),i=1,nbank) close(icsa_history) 850 format(10f8.3) 851 format(5e15.6) ifar=nseed/4*3+1 ifar=nseed+1 ENDIF finished=.false. iter = 0 irecv = 0 isent =0 ifrom= 0 time0i=MPI_WTIME() time1i=time0i time_start_c=time0i if (.not.sync_iter) then time_start=time0i nft00=nft else sync_iter=.false. endif nft00_c=nft nft0i=nft !cccccccccccccccccccccccccccccccccccccc do while (.not. finished) !cccccccccccccccccccccccccccccccccccccc !rc print *,"iter ", iter,' isent=',isent IF (ME.EQ.KING) THEN ! start energy minimization if (isent.eq.0) then ! king-emperor - select seeds & make var & feed input !d print *,'generating new conf',ntrial,MPI_WTIME() call select_is(nseed,ifar,idum) open(icsa_seed,file=csa_seed,status="old") write(icsa_seed,39) & jlee,icycle,nstep,(is(i),bene(is(i)),i=1,nseed) close(icsa_seed) call history_append write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,& ebmin,ebmax,nft,iuse,nbank,ntbank close(icsa_history) call make_var(ntry,idum,iter) !d print *,'new trial generated',ntrial,MPI_WTIME() time2i=MPI_WTIME() write (iout,'(a20,i4,f12.2)') & 'Time for make trial',iter+1,time2i-time1i endif !rc write(iout,*)'1:Calling FEEDIN NTRY',NTRY,' ntrial',ntrial !rc call feedin(ntry,nft) isent=isent+1 if (isent.ge.nodes.or.iter.gt.0) then !t print *,'waiting ',MPI_WTIME() irecv=irecv+1 call recv(0,ifrom,xout,eout,ind,timeout) !t print *,' ',irecv,' received from',ifrom,MPI_WTIME() else ifrom=ifrom+1 endif !t print *,'sending to',ifrom,MPI_WTIME() call send(isent,ifrom,iter) !t print *,isent,' sent ',MPI_WTIME() ! store results ----------------------------------------------- if (isent.ge.nodes.or.iter.gt.0) then nft=nft+ind(3) movernx(irecv)=iabs(ind(5)) call getx(ind,xout,eout,cout,rad,iw_pdb,irecv) if(vdisulf) then nss_out(irecv)=nss do i=1,nss iss_out(i,irecv)=ihpb(i) jss_out(i,irecv)=jhpb(i) enddo endif if(iw_pdb.gt.0) & call write_csa_pdb(xout,eout,nft,irecv,iw_pdb) endif !-------------------------------------------------------------- if (isent.eq.ntry) then time1i=MPI_WTIME() write (iout,'(a18,f12.2,a14,f10.2)') & 'Nonsetup time ',time1i-time_start_c,& ' sec, Eval/s =',(nft-nft00_c)/(time1i-time_start_c) write (iout,'(a14,i4,f12.2,a14,f10.2)') & 'Time for iter ',iter+1,time1i-time0i,& ' sec, Eval/s =',(nft-nft0i)/(time1i-time0i) time0i=time1i nft0i=nft cutdif=cutdif*xctdif if(cutdif.lt.ctdif1) cutdif=ctdif1 if (iter.eq.0) then print *,'UPDATING ',ntry-nodes+1,irecv write(iout,*) 'UPDATING ',ntry-nodes+1 iter=iter+1 !----------------- call update(ntry-nodes+1) ------------------- nstep=nstep+ntry-nseed-(nodes-1) call refresh_bank(ntry-nodes+1) !!bankt call refresh_bankt(ntry-nodes+1) else !----------------- call update(ntry) --------------------------- iter=iter+1 print *,'UPDATING ',ntry,irecv write(iout,*) 'UPDATING ',ntry nstep=nstep+ntry-nseed call refresh_bank(ntry) !!bankt call refresh_bankt(ntry) endif !----------------------------------------------------------------- call write_bank(jlee,nft) !!bankt call write_bankt(jlee,nft) call find_min time1i=MPI_WTIME() write (iout,'(a20,i4,f12.2)') & 'Time for refresh ',iter,time1i-time0i if(ebmin.lt.estop) finished=.true. if(icycle.gt.icmax) then call write_bank1(jlee) do i=1,nbank ibank(i)=2 ibank(i)=1 enddo nbank=nbank+nconf if(nbank.gt.1000) then finished=.true. else !rc goto 333 sync_iter=.true. endif endif if(nstep.gt.nstmax) finished=.true. if(finished.or.sync_iter) then do ij=1,nodes-1 call recv(1,ifrom,xout,eout,ind,timeout) if (timeout) then nstep=nstep+ij-1 print *,'ERROR worker is not responding' write(iout,*) 'ERROR worker is not responding' time1i=MPI_WTIME()-time_start_c print *,'End of cycle, master time for ',iter,' iters ',& time1i,'sec, Eval/s ',(nft-nft00_c)/time1i write (iout,*) 'End of cycle, master time for ',iter,& ' iters ',time1i,' sec' write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i print *,'UPDATING ',ij-1 write(iout,*) 'UPDATING ',ij-1 call flush(iout) call refresh_bank(ij-1) !!bankt call refresh_bankt(ij-1) goto 1002 endif ! print *,'node ',ifrom,' finished ',ij,nft write(iout,*) 'node ',ifrom,' finished ',ij,nft call flush(iout) nft=nft+ind(3) movernx(ij)=iabs(ind(5)) call getx(ind,xout,eout,cout,rad,iw_pdb,ij) if(vdisulf) then nss_out(ij)=nss do i=1,nss iss_out(i,ij)=ihpb(i) jss_out(i,ij)=jhpb(i) enddo endif if(iw_pdb.gt.0) & call write_csa_pdb(xout,eout,nft,ij,iw_pdb) enddo nstep=nstep+nodes-1 !rc print *,'---------bcast finished--------',finished time1i=MPI_WTIME()-time_start_c print *,'End of cycle, master time for ',iter,' iters ',& time1i,'sec, Eval/s ',(nft-nft00_c)/time1i write (iout,*) 'End of cycle, master time for ',iter,& ' iters ',time1i,' sec' write (iout,*) 'Total eval/s ',(nft-nft00_c)/time1i !timeout call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr) !timeout call mpi_bcast(sync_iter,1,mpi_logical,0, !timeout & CG_COMM,ierr) do ij=1,nodes-1 tstart=MPI_WTIME() call mpi_issend(finished,1,mpi_logical,ij,idchar,& CG_COMM,ireq,ierr) call mpi_issend(sync_iter,1,mpi_logical,ij,idchar,& CG_COMM,ireq2,ierr) flag=.false. timeout1=.false. do while(.not. (flag .or. timeout1)) call MPI_TEST(ireq2,flag,muster,ierr) tend1=MPI_WTIME() if(tend1-tstart.gt.60) then print *,'ERROR worker ',ij,' is not responding' write(iout,*) 'ERROR worker ',ij,' is not responding' timeout1=.true. endif enddo if(timeout1) then write(iout,*) 'worker ',ij,' NOT OK ',tend1-tstart timeout=.true. else write(iout,*) 'worker ',ij,' OK ',tend1-tstart endif enddo print *,'UPDATING ',nodes-1,ij write(iout,*) 'UPDATING ',nodes-1 call refresh_bank(nodes-1) !!bankt call refresh_bankt(nodes-1) 1002 continue call write_bank(jlee,nft) !!bankt call write_bankt(jlee,nft) call find_min do i=0,mxmv do j=1,3 nstatnx_tot(i,j)=nstatnx_tot(i,j)+nstatnx(i,j) nstatnx(i,j)=0 enddo enddo write(iout,*)'### Total stats:' do i=0,mxmv if(nstatnx_tot(i,1).ne.0) then if (i.le.9) then write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') & '### N',i,' total=',nstatnx_tot(i,1),& ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',& (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1) else write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') & '###N',i,' total=',nstatnx_tot(i,1),& ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),'%acc',& (nstatnx_tot(i,2)+nstatnx_tot(i,3))*100.0/nstatnx_tot(i,1) endif else if (i.le.9) then write(iout,'(a5,i1,a7,i6,a7,i4,a5,i4,a5,f5.1)') & '### N',i,' total=',nstatnx_tot(i,1),& ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),& ' %acc',0.0 else write(iout,'(a4,i2,a7,i6,a7,i4,a5,i4,a5,f5.1)') & '###N',i,' total=',nstatnx_tot(i,1),& ' close=',nstatnx_tot(i,2),' far=',nstatnx_tot(i,3),& ' %acc',0.0 endif endif enddo endif if(sync_iter) goto 331 39 format(2i3,i7,5(i4,f8.3,1x),19(/,13x,5(i4,f8.3,1x))) 40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4) 43 format(10i8) 44 format('jlee =',i3,':',4f10.1,' E =',f8.3,i7,i10) isent=0 irecv=0 endif ELSE ! soldier - perform energy minimization call minim_jlee print *,'End of minim, proc',me,'time ',MPI_WTIME()-time_start write (iout,*) 'End of minim, proc',me,'time ',& MPI_WTIME()-time_start call flush(iout) !timeout call mpi_bcast(finished,1,mpi_logical,0,CG_COMM,ierr) !timeout call mpi_bcast(sync_iter,1,mpi_logical,0,CG_COMM,ierr) call mpi_recv(finished,1,mpi_logical,0,idchar,& CG_COMM,muster,ierr) call mpi_recv(sync_iter,1,mpi_logical,0,idchar,& CG_COMM,muster,ierr) if(sync_iter) goto 331 ENDIF !cccccccccccccccccccccccccccccccccccccc enddo !cccccccccccccccccccccccccccccccccccccc IF (ME.EQ.KING) THEN call history_append write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,& ebmin,ebmax,nft,iuse,nbank,ntbank write(icsa_history,44) jlee,0.0,0.0,0.0,& 0.0,ebmin,nstep,nft write(icsa_history,*) close(icsa_history) time1i=MPI_WTIME()-time_start print *,'End of RUN, master time ',& time1i,'sec, Eval/s ',(nft-nft00)/time1i write (iout,*) 'End of RUN, master time ',& time1i,' sec' write (iout,*) 'Total eval/s ',(nft-nft00)/time1i if(timeout) then write(iout,*) '!!!! ERROR worker was not responding' write(iout,*) '!!!! cannot finish work normally' write(iout,*) 'Processor0 is calling MPI_ABORT' print *,'!!!! ERROR worker was not responding' print *,'!!!! cannot finish work normally' print *,'Processor0 is calling MPI_ABORT' call flush(iout) call mpi_abort(mpi_comm_world, 111, ierr) endif ENDIF !ccccccccccccccccccccccccccccc 300 continue !ccccccccccccccccccccccccccccc return end subroutine together !----------------------------------------------------------------------------- subroutine feedin(nconf,nft) use MPI_data use geometry_data, only:nvar use io, only:write_csa_pdb ! sends out starting conformations and receives results of energy minimization ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTROL' include 'mpif.h' real(kind=8),dimension(6*nres) :: xin,xout !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout real(kind=8),dimension(2) :: cout integer,dimension(9) :: ind integer,dimension(12) :: info integer,dimension(mpi_status_size) :: muster ! include 'COMMON.SETUP' real(kind=8),parameter :: rad=1.745329252d-2 integer :: j,nconf,nft,mm,n,ierror,ierrcode,ierr,iw_pdb,& man print *,'FEEDIN: NCONF=',nconf mm=0 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (nconf .lt. nodes-1) then write (*,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',& nconf,nodes-1 write (iout,*) 'FATAL ERROR in FEEDIN, nconf < nodes -1',& nconf,nodes-1 call mpi_abort(mpi_comm_world,ierror,ierrcode) endif do n=1,nconf ! pull out external and internal variables for next start call putx(xin,n,rad) ! write (iout,*) 'XIN from FEEDIN N=',n ! write(iout,'(8f10.4)') (xin(j),j=1,nvar) mm=mm+1 if (mm.lt.nodes) then ! feed task to soldier ! print *, ' sending input for start # ',n info(1)=n info(2)=-1 info(3)=0 info(4)=0 info(5)=0 info(6)=0 call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,& ierr) call mpi_send(xin,nvar,mpi_double_precision,mm,& idreal,CG_COMM,ierr) else ! find an available soldier call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,& CG_COMM,muster,ierr) ! print *, ' receiving output from start # ',ind(1) man=muster(mpi_source) ! receive final energies and variables nft=nft+ind(3) call mpi_recv(eout,1,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) ! print *,eout #ifdef CO_BIAS call mpi_recv(co,1,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co #endif call mpi_recv(xout,nvar,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) ! print *,nvar , ierr ! feed next task to soldier ! print *, ' sending input for start # ',n info(1)=n info(2)=-1 info(3)=0 info(4)=0 info(5)=0 info(6)=0 info(7)=0 info(8)=0 info(9)=0 call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,& ierr) call mpi_send(xin,nvar,mpi_double_precision,man,& idreal,CG_COMM,ierr) ! retrieve latest results call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1)) if(iw_pdb.gt.0) & call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb) endif enddo !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! no more input !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do j=1,nodes-1 ! wait for a soldier call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint,& CG_COMM,muster,ierr) !rc if (ierr.ne.0) go to 30 ! print *, ' receiving output from start # ',ind(1) man=muster(mpi_source) ! receive final energies and variables nft=nft+ind(3) call mpi_recv(eout,1,& mpi_double_precision,man,idreal,& CG_COMM,muster,ierr) ! print *,eout #ifdef CO_BIAS call mpi_recv(co,1,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co #endif !rc if (ierr.ne.0) go to 30 call mpi_recv(xout,nvar,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) ! print *,nvar , ierr !rc if (ierr.ne.0) go to 30 ! halt soldier info(1)=0 info(2)=-1 info(3)=0 info(4)=0 info(5)=0 info(6)=0 call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,& ierr) ! retrieve results call getx(ind,xout,eout,cout,rad,iw_pdb,ind(1)) if(iw_pdb.gt.0) & call write_csa_pdb(xout,eout,nft,ind(1),iw_pdb) enddo !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc return 10 print *, ' dispatching error' call mpi_abort(mpi_comm_world,ierror,ierrcode) return 20 print *, ' communication error' call mpi_abort(mpi_comm_world,ierror,ierrcode) return 30 print *, ' receiving error' call mpi_abort(mpi_comm_world,ierror,ierrcode) return end subroutine feedin !----------------------------------------------------------------------------- subroutine getx(ind,xout,eout,cout,rad,iw_pdb,k) use geometry_data use energy_data use compare, only: contact_fract use MPI_data include 'mpif.h' ! receives and stores data from soldiers ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.CONTACTS' integer,dimension(9) :: ind real(kind=8),dimension(6*nres) :: xout !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout !jlee real(kind=8) :: przes(3),obr(3,3),cout(2) logical :: non_conv integer :: iw_pdb,k,j,ierror,ierrcode real(kind=8) :: rad,co !jlee iw_pdb=2 if (k.gt.mxio .or. k.lt.1) then write (iout,*) & 'ERROR - dimensions of ANGMIN have been exceeded K=',k call mpi_abort(mpi_comm_world,ierror,ierrcode) endif ! store ind() do j=1,9 indb(k,j)=ind(j) enddo ! store energies etot(k)=eout(1) ! retrieve dihedral angles etc call var_to_geom(nvar,xout) do j=2,nres-1 dihang(1,j,1,k)=theta(j+1) dihang(2,j,1,k)=phi(j+2) dihang(3,j,1,k)=alph(j) dihang(4,j,1,k)=omeg(j) enddo dihang(2,nres-1,1,k)=0.0d0 !jlee if(iref.eq.0) then iw_pdb=1 !d write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a4,i3,i4)') !d & ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' mv ', !d & ind(5),ind(4) return endif call chainbuild ! call dihang_to_c(dihang(1,1,1,k)) ! call fitsq(rms,c(1,1),crefjlee(1,1),nres,przes,obr,non_conv) ! call fitsq(rms,c(1,2),crefjlee(1,2),nres-1,przes,obr,non_conv) ! call fitsq(rms,c(1,nstart_seq),crefjlee(1,nstart_sup), ! & nsup,przes,obr,non_conv) ! rmsn(k)=dsqrt(rms) call rmsd_csa(rmsn(k)) call contact(.false.,ncont,icont,co) pncn(k)=contact_fract(ncont,ncont_ref,icont,icont_ref) !d write(iout,'(i3,a3,i4,i5,a6,1pe12.4,a5 !d & ,0pf5.2,a5,f5.1,a,f6.3,a4,i3,i4)') !d & ind(2),' e ',ind(3),ind(1),' etot ',etot(k),' rms ', !d & rmsn(k),' %NC ',pncn(k)*100,' cont.order',co,' mv ', !d & ind(5),ind(4) if (rmsn(k).gt.rmscut.or.pncn(k).lt.pnccut) iw_pdb=0 return end subroutine getx !----------------------------------------------------------------------------- subroutine putx(xin,n,rad) use geometry_data ! gets starting variables ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' integer :: n,m,j real(kind=8),dimension(6*nres) :: xin !(maxvar) (maxvar=6*maxres) real(kind=8) :: rad ! pull out starting values for variables ! write (iout,*)'PUTX: N=',n do m=1,numch ! write (iout,'(8f10.4)') (dihang_in(1,j,m,n),j=2,nres-1) ! write (iout,'(8f10.4)') (dihang_in(2,j,m,n),j=2,nres-1) ! write (iout,'(8f10.4)') (dihang_in(3,j,m,n),j=2,nres-1) ! write (iout,'(8f10.4)') (dihang_in(4,j,m,n),j=2,nres-1) do j=2,nres-1 theta(j+1)=dihang_in(1,j,m,n) phi(j+2)=dihang_in(2,j,m,n) alph(j)=dihang_in(3,j,m,n) omeg(j)=dihang_in(4,j,m,n) enddo enddo ! set up array of variables call geom_to_var(nvar,xin) ! write (iout,*) 'xin in PUTX N=',n ! call intout ! write (iout,'(8f10.4)') (xin(i),i=1,nvar) return end subroutine putx !----------------------------------------------------------------------------- subroutine putx2(xin,iff,n) use geometry_data ! gets starting variables ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' integer :: n,m,j,i real(kind=8),dimension(6*nres) :: xin !(maxvar) (maxvar=6*maxres) integer,dimension(nres) :: iff !(maxres) ! pull out starting values for variables do m=1,numch do j=2,nres-1 theta(j+1)=dihang_in2(1,j,m,n) phi(j+2)=dihang_in2(2,j,m,n) alph(j)=dihang_in2(3,j,m,n) omeg(j)=dihang_in2(4,j,m,n) enddo enddo ! set up array of variables call geom_to_var(nvar,xin) do i=1,nres iff(i)=iff_in(i,n) enddo return end subroutine putx2 !----------------------------------------------------------------------------- subroutine prune_bank(p_cut) use MPI_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.TIME1' ! include 'COMMON.SETUP' integer :: k,j,i,m,ip,nprune real(kind=8) :: p_cut,diff,ddmin !--------------------------- ! This subroutine prunes bank conformations using p_cut !--------------------------- nprune=0 nprune=nprune+1 m=1 do k=1,numch do j=2,nres-1 do i=1,4 dihang(i,j,k,nprune)=bvar(i,j,k,m) enddo enddo enddo bene(nprune)=bene(m) brmsn(nprune)=brmsn(m) bpncn(nprune)=bpncn(m) do m=2,nbank ddmin=9.d190 do ip=1,nprune call get_diff12(dihang(1,1,1,ip),bvar(1,1,1,m),diff) if(diff.lt.p_cut) goto 100 if(diff.lt.ddmin) ddmin=diff enddo nprune=nprune+1 do k=1,numch do j=2,nres-1 do i=1,4 dihang(i,j,k,nprune)=bvar(i,j,k,m) enddo enddo enddo bene(nprune)=bene(m) brmsn(nprune)=brmsn(m) bpncn(nprune)=bpncn(m) 100 continue write (iout,*) 'Pruning :',m,nprune,p_cut,ddmin enddo nbank=nprune print *, 'Pruning :',m,nprune,p_cut call write_bank(0,0) return end subroutine prune_bank !----------------------------------------------------------------------------- subroutine reminimize(jlee) use MPI_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' include 'mpif.h' ! include 'COMMON.CSA' ! include 'COMMON.BANK' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.TIME1' ! include 'COMMON.SETUP' integer :: i,j,k,jlee,index,nft,ntry !--------------------------- ! This subroutine re-minimizes bank conformations: !--------------------------- ntry=nbank call find_max call find_min if (me.eq.king) then open(icsa_history,file=csa_history,status="old") write(icsa_history,*) "Re-minimization",nodes,"nodes" write(icsa_history,851) (bene(i),i=1,nbank) write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,& ebmin,ebmax,nft,iuse,nbank,ntbank close(icsa_history) do index=1,ntry do k=1,numch do j=2,nres-1 do i=1,4 dihang_in(i,j,k,index)=bvar(i,j,k,index) enddo enddo enddo enddo nft=0 call feedin(ntry,nft) else call minim_jlee endif call find_max call find_min if (me.eq.king) then do i=1,ntry call replace_bvar(i,i) enddo open(icsa_history,file=csa_history,status="old") write(icsa_history,40) jlee,icycle,nstep,cutdif,ibmin,ibmax,& ebmin,ebmax,nft,iuse,nbank,ntbank write(icsa_history,851) (bene(i),i=1,nbank) close(icsa_history) call write_bank_reminimized(jlee,nft) endif 40 format(2i2,i8,f8.1,2i4,2(1pe14.5),i10,3i4) 851 format(5e15.6) 850 format(5e15.10) ! 850 format(10f8.3) return end subroutine reminimize !----------------------------------------------------------------------------- subroutine send(n,mm,it) use MPI_data use geometry_data, only: nvar use control_data, only: vdisulf ! sends out starting conformation for minimization ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTROL' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' include 'mpif.h' real(kind=8),dimension(6*nres) :: xin,xout,xin2 !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout real(kind=8),dimension(2) :: cout integer,dimension(9) :: ind integer,dimension(nres) :: iff !(maxres) integer,dimension(12) :: info integer,dimension(mpi_status_size) :: muster ! include 'COMMON.SETUP' real(kind=8),parameter :: rad=1.745329252d-2 integer :: n,mm,it,ierr if (isend2(n).eq.0) then ! pull out external and internal variables for next start call putx(xin,n,rad) info(1)=n info(2)=it info(3)=movenx(n) info(4)=nss_in(n) info(5)=parent(1,n) info(6)=parent(2,n) if (movenx(n).eq.14.or.movenx(n).eq.17) then info(7)=idata(1,n) info(8)=idata(2,n) else if (movenx(n).eq.16) then info(7)=idata(1,n) info(8)=idata(2,n) info(10)=idata(3,n) info(11)=idata(4,n) info(12)=idata(5,n) else info(7)=0 info(8)=0 info(10)=0 info(11)=0 info(12)=0 endif if (movenx(n).eq.15) then info(9)=parent(3,n) else info(9)=0 endif call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,& ierr) call mpi_send(xin,nvar,mpi_double_precision,mm,& idreal,CG_COMM,ierr) else ! distfit & minimization for n7 move info(1)=-n info(2)=it info(3)=movenx(n) info(4)=nss_in(n) info(5)=parent(1,n) info(6)=parent(2,n) info(7)=0 info(8)=0 info(9)=0 call mpi_send(info,12,mpi_integer,mm,idint,CG_COMM,& ierr) call putx2(xin,iff,isend2(n)) call mpi_send(xin,nvar,mpi_double_precision,mm,& idreal,CG_COMM,ierr) call mpi_send(iff,nres,mpi_integer,mm,& idint,CG_COMM,ierr) call putx(xin2,n,rad) call mpi_send(xin2,nvar,mpi_double_precision,mm,& idreal,CG_COMM,ierr) endif if (vdisulf.and.nss_in(n).ne.0) then call mpi_send(iss_in(1,n),nss_in(n),mpi_integer,mm,& idint,CG_COMM,ierr) call mpi_send(jss_in(1,n),nss_in(n),mpi_integer,mm,& idint,CG_COMM,ierr) endif return end subroutine send !----------------------------------------------------------------------------- subroutine recv(ihalt,man,xout,eout,ind,tout) use MPI_data use energy_data use geometry_data, only: nvar use control_data, only: vdisulf ! receives results of energy minimization ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTROL' ! include 'COMMON.SBRIDGE' ! include 'COMMON.BANK' ! include 'COMMON.CHAIN' include 'mpif.h' real(kind=8),dimension(6*nres) :: xin,xout !(maxvar) (maxvar=6*maxres) real(kind=8),dimension(mxch*(mxch+1)/2+1) :: eout real(kind=8),dimension(2) :: cout integer,dimension(9) :: ind integer,dimension(12) :: info integer,dimension(mpi_status_size) :: muster ! include 'COMMON.SETUP' logical :: tout,flag real(kind=8) :: tstart,tend1 real(kind=8),parameter :: twait=600.0d0 integer :: ihalt,man,ierr ! find an available soldier tout=.false. flag=.false. tstart=MPI_WTIME() do while(.not. (flag .or. tout)) call MPI_IPROBE(mpi_any_source,idint,CG_COMM,flag, & muster,ierr) tend1=MPI_WTIME() if(tend1-tstart.gt.twait .and. ihalt.eq.1) tout=.true. !_error if(tend1-tstart.gt.twait) tout=.true. enddo if (tout) then write(iout,*) 'ERROR = timeout for recv ',tend1-tstart call flush(iout) return endif man=muster(mpi_source) !timeout call mpi_recv(ind,9,mpi_integer,mpi_any_source,idint, !timeout * CG_COMM,muster,ierr) ! print *, ' receiving output from start # ',ind(1) !t print *,'receiving ',MPI_WTIME() !timeout man=muster(mpi_source) call mpi_recv(ind,9,mpi_integer,man,idint,& CG_COMM,muster,ierr) !timeout ! receive final energies and variables call mpi_recv(eout,1,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) ! print *,eout #ifdef CO_BIAS call mpi_recv(co,1,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) write (iout,'(a15,f3.2,$)') ' BIAS by contact order*100 ',co #endif call mpi_recv(xout,nvar,mpi_double_precision,& man,idreal,CG_COMM,muster,ierr) ! print *,nvar , ierr if(vdisulf) nss=ind(6) if(vdisulf.and.nss.ne.0) then call mpi_recv(ihpb,nss,mpi_integer,& man,idint,CG_COMM,muster,ierr) call mpi_recv(jhpb,nss,mpi_integer,& man,idint,CG_COMM,muster,ierr) endif ! halt soldier if(ihalt.eq.1) then ! print *,'sending halt to ',man write(iout,*) 'sending halt to ',man info(1)=0 call mpi_send(info,12,mpi_integer,man,idint,CG_COMM,ierr) endif return end subroutine recv !----------------------------------------------------------------------------- subroutine history_append ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' #if defined(AIX) || defined(PGI) open(icsa_history,file=csa_history,position="append") #else open(icsa_history,file=csa_history,access="append") #endif return end subroutine history_append !----------------------------------------------------------------------------- subroutine alloc_CSA_arrays use energy_data, only: ns mxgr=2*nres if(.not.allocated(bfrag)) allocate(bfrag(4,nres/3)) ! commom.bank ! common/varin/ !el allocate(dihang_in(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio) allocate(dihang_in(mxang,nres,mxch,5000)) !(mxang,maxres,mxch,mxio) allocate(nss_in(mxio)) !(mxio) allocate(iss_in(ns,mxio),jss_in(ns,mxio)) !(maxss,mxio) ! common/minvar/ allocate(dihang(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio) allocate(rmsn(mxio),pncn(mxio)) !(mxio) allocate(etot(mxio)) !(mxio) allocate(nss_out(mxio)) !(mxio) allocate(iss_out(ns,mxio),jss_out(ns,mxio)) !(maxss,mxio) ! common/bank/ allocate(rvar(mxang,nres,mxch,mxio),bvar(mxang,nres,mxch,mxio)) !(mxang,maxres,mxch,mxio) allocate(bene(mxio),rene(mxio),brmsn(mxio),rrmsn(mxio)) allocate(bpncn(mxio),rpncn(mxio)) !(mxio) allocate(ibank(mxio),is(mxio),jbank(mxio)) !(mxio) allocate(dij(mxio,mxio)) !(mxio,mxio) ! common/bank_disulfid/ allocate(bvar_nss(mxio),bvar_ns(mxio)) !(mxio) allocate(bvar_s(ns,mxio)) !(maxss,mxio) allocate(bvar_ss(2,ns,mxio)) !(2,maxss,mxio) ! common/mvstat/ allocate(movenx(mxio),movernx(mxio)) !(mxio) allocate(nstatnx(0:mxmv,3),nstatnx_tot(0:mxmv,3)) !(0:mxmv,3) allocate(indb(mxio,9)) !(mxio,9) allocate(parent(3,mxio)) !(3,mxio) ! common/send2/ allocate(isend2(mxio)) !(mxio) allocate(iff_in(nres,mxio2)) !(maxres,mxio2) allocate(dihang_in2(mxang,nres,mxch,mxio2)) !(mxang,maxres,mxch,mxio2) allocate(idata(5,mxio)) !(5,mxio) ! common.csa ! common/alphaa/ allocate(ngroup(mxgr)) !(mxgr) allocate(igroup(3,mxang,mxgr)) !(3,mxang,mxgr) ! common.distfit ! COMMON /frag/ allocate(bvar_frag(mxio,6)) !(mxio,6) allocate(hvar_frag(mxio,3),lvar_frag(mxio,3),svar_frag(mxio,3)) !(mxio,3) allocate(avar_frag(mxio,5)) !(mxio,5) ! commom.hairpin ! common /spinka/ allocate(nharp_seed(nseed),nharp_use(nseed)) !(max_seed) allocate(iharp_seed(4,nres/3,nseed)) !(4,maxres/3,max_seed) allocate(iharp_use(0:4,nres/3,nseed)) !(0:4,maxres/3,max_seed) return end subroutine alloc_CSA_arrays !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module csa