--- /dev/null
+ subroutine mcm_setup
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MCM'
+ include 'COMMON.CONTROL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+C
+C Set up variables used in MC/MCM.
+C
+ write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
+ write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,
+ & ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
+ write (iout,'(4(a,f8.1)/2(a,i3))')
+ & 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,
+ & ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
+ if (nwindow.gt.0) then
+ write (iout,'(a)') 'Perturbation windows:'
+ do i=1,nwindow
+ i1=winstart(i)
+ i2=winend(i)
+ it1=itype(i1)
+ it2=itype(i2)
+ write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,
+ & ' length',winlen(i)
+ enddo
+ endif
+C Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
+ RBol=1.9858D-3
+C Number of "end bonds".
+ koniecl=0
+c koniecl=nphi
+ print *,'koniecl=',koniecl
+ write (iout,'(a)') 'Probabilities of move types:'
+ write (*,'(a)') 'Probabilities of move types:'
+ do i=1,MaxMoveType
+ write (iout,'(a,f10.5)') MovTypID(i),
+ & sumpro_type(i)-sumpro_type(i-1)
+ write (*,'(a,f10.5)') MovTypID(i),
+ & sumpro_type(i)-sumpro_type(i-1)
+ enddo
+ write (iout,*)
+C Maximum length of N-bond segment to be moved
+c nbm=nres-1-(2*koniecl-1)
+ if (nwindow.gt.0) then
+ maxwinlen=winlen(1)
+ do i=2,nwindow
+ if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
+ enddo
+ nbm=min0(maxwinlen,6)
+ write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
+ else
+ nbm=min0(6,nres-2)
+ endif
+ sumpro_bond(0)=0.0D0
+ sumpro_bond(1)=0.0D0
+ do i=2,nbm
+ sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
+ enddo
+ write (iout,'(a)') 'The SumPro_Bond array:'
+ write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
+ write (*,'(a)') 'The SumPro_Bond array:'
+ write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
+C Maximum number of side chains moved simultaneously
+c print *,'nnt=',nnt,' nct=',nct
+ ngly=0
+ do i=nnt,nct
+ if (itype(i).eq.10) ngly=ngly+1
+ enddo
+ mmm=nct-nnt-ngly+1
+ if (mmm.gt.0) then
+ MaxSideMove=min0((nct-nnt+1)/2,mmm)
+ endif
+c print *,'MaxSideMove=',MaxSideMove
+C Max. number of generated confs (not used at present).
+ maxgen=10000
+C Set initial temperature
+ Tcur=Tmin
+ betbol=1.0D0/(Rbol*Tcur)
+ write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,
+ & ' BetBol:',betbol
+ write (iout,*) 'RanFract=',ranfract
+ return
+ end
+c------------------------------------------------------------------------------
+#ifndef MPI
+ subroutine do_mcm(i_orig)
+C Monte-Carlo-with-Minimization calculations - serial code.
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.CHAIN'
+ include 'COMMON.MCM'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.CACHE'
+crc include 'COMMON.DEFORM'
+crc include 'COMMON.DEFORM1'
+ include 'COMMON.NAMES'
+ logical accepted,over,ovrtim,error,lprint,not_done,my_conf,
+ & enelower,non_conv
+ integer MoveType,nbond,conf_comp
+ integer ifeed(max_cache)
+ double precision varia(maxvar),varold(maxvar),elowest,eold,
+ & przes(3),obr(3,3)
+ double precision energia(0:n_ene)
+ double precision coord1(maxres,3)
+
+C---------------------------------------------------------------------------
+C Initialize counters.
+C---------------------------------------------------------------------------
+C Total number of generated confs.
+ ngen=0
+C Total number of moves. In general this won't be equal to the number of
+C attempted moves, because we may want to reject some "bad" confs just by
+C overlap check.
+ nmove=0
+C Total number of temperature jumps.
+ ntherm=0
+C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
+C motions.
+ ncache=0
+ do i=1,nres
+ nbond_move(i)=0
+ enddo
+C Initialize total and accepted number of moves of various kind.
+ do i=0,MaxMoveType
+ moves(i)=0
+ moves_acc(i)=0
+ enddo
+C Total number of energy evaluations.
+ neneval=0
+ nfun=0
+ nsave=0
+
+ write (iout,*) 'RanFract=',RanFract
+
+ WhatsUp=0
+ Kwita=0
+
+c----------------------------------------------------------------------------
+C Compute and print initial energies.
+c----------------------------------------------------------------------------
+ call intout
+ write (iout,'(/80(1h*)/a)') 'Initial energies:'
+ call chainbuild
+ nf=0
+
+ call etotal(energia(0))
+ etot = energia(0)
+C Minimize the energy of the first conformation.
+ if (minim) then
+ call geom_to_var(nvar,varia)
+! write (iout,*) 'The VARIA array'
+! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
+ call minimize(etot,varia,iretcode,nfun)
+ call var_to_geom(nvar,varia)
+ call chainbuild
+ write (iout,*) 'etot from MINIMIZE:',etot
+! write (iout,*) 'Tha VARIA array'
+! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
+
+ call etotal(energia(0))
+ etot=energia(0)
+ call enerprint(energia(0))
+ endif
+ if (refstr) then
+ call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
+ & obr,non_conv)
+ rms=dsqrt(rms)
+ call contact(.false.,ncont,icont,co)
+ frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
+ write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
+ & 'RMS deviation from the reference structure:',rms,
+ & ' % of native contacts:',frac*100,' contact order:',co
+ if (print_stat)
+ & write (istat,'(i5,17(1pe14.5))') 0,
+ & (energia(print_order(i)),i=1,nprint_ene),
+ & etot,rms,frac,co
+ else
+ if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,
+ & (energia(print_order(i)),i=1,nprint_ene),etot
+ endif
+ if (print_stat) close(istat)
+ neneval=neneval+nfun+1
+ write (iout,'(/80(1h*)/20x,a/80(1h*))')
+ & 'Enter Monte Carlo procedure.'
+ if (print_int) then
+ close(igeom)
+ call briefout(0,etot)
+ endif
+ eold=etot
+ do i=1,nvar
+ varold(i)=varia(i)
+ enddo
+ elowest=etot
+ call zapis(varia,etot)
+ nacc=0 ! total # of accepted confs of the current processor.
+ nacc_tot=0 ! total # of accepted confs of all processors.
+
+ not_done = (iretcode.ne.11)
+
+C----------------------------------------------------------------------------
+C Main loop.
+c----------------------------------------------------------------------------
+ it=0
+ nout=0
+ do while (not_done)
+ it=it+1
+ write (iout,'(80(1h*)/20x,a,i7)')
+ & 'Beginning iteration #',it
+C Initialize local counter.
+ ntrial=0 ! # of generated non-overlapping confs.
+ accepted=.false.
+ do while (.not. accepted)
+
+C Retrieve the angles of previously accepted conformation
+ noverlap=0 ! # of overlapping confs.
+ do j=1,nvar
+ varia(j)=varold(j)
+ enddo
+ call var_to_geom(nvar,varia)
+C Rebuild the chain.
+ call chainbuild
+C Heat up the system, if necessary.
+ call heat(over)
+C If temperature cannot be further increased, stop.
+ if (over) goto 20
+ MoveType=0
+ nbond=0
+ lprint=.true.
+cd write (iout,'(a)') 'Old variables:'
+cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+C Decide whether to generate a random conformation or perturb the old one
+ RandOrPert=ran_number(0.0D0,1.0D0)
+ if (RandOrPert.gt.RanFract) then
+ if (print_mc.gt.0)
+ & write (iout,'(a)') 'Perturbation-generated conformation.'
+ call perturb(error,lprint,MoveType,nbond,1.0D0)
+ if (error) goto 20
+ if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
+ write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
+ & MoveType,' returned from PERTURB.'
+ goto 20
+ endif
+ call chainbuild
+ else
+ MoveType=0
+ moves(0)=moves(0)+1
+ nstart_grow=iran_num(3,nres)
+ if (print_mc.gt.0)
+ & write (iout,'(2a,i3)') 'Random-generated conformation',
+ & ' - chain regrown from residue',nstart_grow
+ call gen_rand_conf(nstart_grow,*30)
+ endif
+ call geom_to_var(nvar,varia)
+cd write (iout,'(a)') 'New variables:'
+cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+ ngen=ngen+1
+
+ call etotal(energia(0))
+ etot=energia(0)
+c call enerprint(energia(0))
+c write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
+ if (etot-elowest.gt.overlap_cut) then
+ if(iprint.gt.1.or.etot.lt.1d20)
+ & write (iout,'(a,1pe14.5)')
+ & 'Overlap detected in the current conf.; energy is',etot
+ neneval=neneval+1
+ accepted=.false.
+ noverlap=noverlap+1
+ if (noverlap.gt.maxoverlap) then
+ write (iout,'(a)') 'Too many overlapping confs.'
+ goto 20
+ endif
+ else
+ if (minim) then
+ call minimize(etot,varia,iretcode,nfun)
+cd write (iout,*) 'etot from MINIMIZE:',etot
+cd write (iout,'(a)') 'Variables after minimization:'
+cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+
+ call etotal(energia(0))
+ etot = energia(0)
+ neneval=neneval+nfun+2
+ endif
+c call enerprint(energia(0))
+ write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,
+ & ' energy:',etot
+C--------------------------------------------------------------------------
+C... Do Metropolis test
+C--------------------------------------------------------------------------
+ accepted=.false.
+ my_conf=.false.
+
+ if (WhatsUp.eq.0 .and. Kwita.eq.0) then
+ call metropolis(nvar,varia,varold,etot,eold,accepted,
+ & my_conf,EneLower)
+ endif
+ write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
+ if (accepted) then
+
+ nacc=nacc+1
+ nacc_tot=nacc_tot+1
+ if (elowest.gt.etot) elowest=etot
+ moves_acc(MoveType)=moves_acc(MoveType)+1
+ if (MoveType.eq.1) then
+ nbond_acc(nbond)=nbond_acc(nbond)+1
+ endif
+C Check against conformation repetitions.
+ irepet=conf_comp(varia,etot)
+ if (print_stat) then
+#if defined(AIX) || defined(PGI)
+ open (istat,file=statname,position='append')
+#else
+ open (istat,file=statname,access='append')
+#endif
+ endif
+ call statprint(nacc,nfun,iretcode,etot,elowest)
+ if (refstr) then
+ call var_to_geom(nvar,varia)
+ call chainbuild
+ call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
+ & nsup,przes,obr,non_conv)
+ rms=dsqrt(rms)
+ call contact(.false.,ncont,icont,co)
+ frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
+ write (iout,'(a,f8.3,a,f8.3)')
+ & 'RMS deviation from the reference structure:',rms,
+ & ' % of native contacts:',frac*100,' contact order',co
+ endif ! refstr
+ if (My_Conf) then
+ nout=nout+1
+ write (iout,*) 'Writing new conformation',nout
+ if (refstr) then
+ write (istat,'(i5,16(1pe14.5))') nout,
+ & (energia(print_order(i)),i=1,nprint_ene),
+ & etot,rms,frac
+ else
+ if (print_stat)
+ & write (istat,'(i5,17(1pe14.5))') nout,
+ & (energia(print_order(i)),i=1,nprint_ene),etot
+ endif ! refstr
+ if (print_stat) close(istat)
+C Print internal coordinates.
+ if (print_int) call briefout(nout,etot)
+C Accumulate the newly accepted conf in the coord1 array, if it is different
+C from all confs that are already there.
+ call compare_s1(n_thr,max_thread2,etot,varia,ii,
+ & enetb1,coord1,rms_deform,.true.,iprint)
+ write (iout,*) 'After compare_ss: n_thr=',n_thr
+ if (ii.eq.1 .or. ii.eq.3) then
+ write (iout,'(8f10.4)')
+ & (rad2deg*coord1(i,n_thr),i=1,nvar)
+ endif
+ else
+ write (iout,*) 'Conformation from cache, not written.'
+ endif ! My_Conf
+
+ if (nrepm.gt.maxrepm) then
+ write (iout,'(a)') 'Too many conformation repetitions.'
+ goto 20
+ endif
+C Store the accepted conf. and its energy.
+ eold=etot
+ do i=1,nvar
+ varold(i)=varia(i)
+ enddo
+ if (irepet.eq.0) call zapis(varia,etot)
+C Lower the temperature, if necessary.
+ call cool
+
+ else
+
+ ntrial=ntrial+1
+ endif ! accepted
+ endif ! overlap
+
+ 30 continue
+ enddo ! accepted
+C Check for time limit.
+ if (ovrtim()) WhatsUp=-1
+ not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
+ & .and. (Kwita.eq.0)
+
+ enddo ! not_done
+ goto 21
+ 20 WhatsUp=-3
+
+ 21 continue
+ runtime=tcpu()
+ write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
+ call statprint(nacc,nfun,iretcode,etot,elowest)
+ write (iout,'(a)')
+ & 'Statistics of multiple-bond motions. Total motions:'
+ write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
+ write (iout,'(a)') 'Accepted motions:'
+ write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
+ if (it.ge.maxacc)
+ &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
+
+ return
+ end
+#endif
+#ifdef MPI
+c------------------------------------------------------------------------------
+ subroutine do_mcm(i_orig)
+C Monte-Carlo-with-Minimization calculations - parallel code.
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'mpif.h'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.CHAIN'
+ include 'COMMON.MCM'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.INFO'
+ include 'COMMON.CACHE'
+crc include 'COMMON.DEFORM'
+crc include 'COMMON.DEFORM1'
+crc include 'COMMON.DEFORM2'
+ include 'COMMON.MINIM'
+ include 'COMMON.NAMES'
+ logical accepted,over,ovrtim,error,lprint,not_done,similar,
+ & enelower,non_conv,flag,finish
+ integer MoveType,nbond,conf_comp
+ double precision varia(maxvar),varold(maxvar),elowest,eold,
+ & x1(maxvar), varold1(maxvar), przes(3),obr(3,3)
+ integer iparentx(max_threadss2)
+ integer iparentx1(max_threadss2)
+ integer imtasks(150),imtasks_n
+ double precision energia(0:n_ene)
+
+ print *,'Master entered DO_MCM'
+ nodenum = nprocs
+
+ finish=.false.
+ imtasks_n=0
+ do i=1,nodenum-1
+ imtasks(i)=0
+ enddo
+C---------------------------------------------------------------------------
+C Initialize counters.
+C---------------------------------------------------------------------------
+C Total number of generated confs.
+ ngen=0
+C Total number of moves. In general this won`t be equal to the number of
+C attempted moves, because we may want to reject some "bad" confs just by
+C overlap check.
+ nmove=0
+C Total number of temperature jumps.
+ ntherm=0
+C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
+C motions.
+ ncache=0
+ do i=1,nres
+ nbond_move(i)=0
+ enddo
+C Initialize total and accepted number of moves of various kind.
+ do i=0,MaxMoveType
+ moves(i)=0
+ moves_acc(i)=0
+ enddo
+C Total number of energy evaluations.
+ neneval=0
+ nfun=0
+ nsave=0
+c write (iout,*) 'RanFract=',RanFract
+ WhatsUp=0
+ Kwita=0
+c----------------------------------------------------------------------------
+C Compute and print initial energies.
+c----------------------------------------------------------------------------
+ call intout
+ write (iout,'(/80(1h*)/a)') 'Initial energies:'
+ call chainbuild
+ nf=0
+ call etotal(energia(0))
+ etot = energia(0)
+ call enerprint(energia(0))
+C Request energy computation from slave processors.
+ call geom_to_var(nvar,varia)
+! write (iout,*) 'The VARIA array'
+! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
+ call minimize(etot,varia,iretcode,nfun)
+ call var_to_geom(nvar,varia)
+ call chainbuild
+ write (iout,*) 'etot from MINIMIZE:',etot
+! write (iout,*) 'Tha VARIA array'
+! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
+ neneval=0
+ eneglobal=1.0d99
+ if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))')
+ & 'Enter Monte Carlo procedure.'
+ if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
+ eold=etot
+ do i=1,nvar
+ varold(i)=varia(i)
+ enddo
+ elowest=etot
+ call zapis(varia,etot)
+c diagnostics
+ call var_to_geom(nvar,varia)
+ call chainbuild
+ call etotal(energia(0))
+ if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
+c end diagnostics
+ nacc=0 ! total # of accepted confs of the current processor.
+ nacc_tot=0 ! total # of accepted confs of all processors.
+ not_done=.true.
+C----------------------------------------------------------------------------
+C Main loop.
+c----------------------------------------------------------------------------
+ it=0
+ nout=0
+ LOOP1:do while (not_done)
+ it=it+1
+ if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
+ & 'Beginning iteration #',it
+C Initialize local counter.
+ ntrial=0 ! # of generated non-overlapping confs.
+ noverlap=0 ! # of overlapping confs.
+ accepted=.false.
+ LOOP2:do while (.not. accepted)
+
+ LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
+ do i=1,nodenum-1
+ if(imtasks(i).eq.0) then
+ is=i
+ exit
+ endif
+ enddo
+C Retrieve the angles of previously accepted conformation
+ do j=1,nvar
+ varia(j)=varold(j)
+ enddo
+ call var_to_geom(nvar,varia)
+C Rebuild the chain.
+ call chainbuild
+C Heat up the system, if necessary.
+ call heat(over)
+C If temperature cannot be further increased, stop.
+ if (over) then
+ finish=.true.
+ endif
+ MoveType=0
+ nbond=0
+c write (iout,'(a)') 'Old variables:'
+c write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+C Decide whether to generate a random conformation or perturb the old one
+ RandOrPert=ran_number(0.0D0,1.0D0)
+ if (RandOrPert.gt.RanFract) then
+ if (print_mc.gt.0)
+ & write (iout,'(a)') 'Perturbation-generated conformation.'
+ call perturb(error,lprint,MoveType,nbond,1.0D0)
+c print *,'after perturb',error,finish
+ if (error) finish = .true.
+ if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
+ write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
+ & MoveType,' returned from PERTURB.'
+ finish=.true.
+ write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',
+ & MoveType,' returned from PERTURB.'
+ endif
+ call chainbuild
+ else
+ MoveType=0
+ moves(0)=moves(0)+1
+ nstart_grow=iran_num(3,nres)
+ if (print_mc.gt.0)
+ & write (iout,'(2a,i3)') 'Random-generated conformation',
+ & ' - chain regrown from residue',nstart_grow
+ call gen_rand_conf(nstart_grow,*30)
+ endif
+ call geom_to_var(nvar,varia)
+ ngen=ngen+1
+c print *,'finish=',finish
+ if (etot-elowest.gt.overlap_cut) then
+ if (print_mc.gt.1) write (iout,'(a,1pe14.5)')
+ & 'Overlap detected in the current conf.; energy is',etot
+ if(iprint.gt.1.or.etot.lt.1d19) print *,
+ & 'Overlap detected in the current conf.; energy is',etot
+ neneval=neneval+1
+ accepted=.false.
+ noverlap=noverlap+1
+ if (noverlap.gt.maxoverlap) then
+ write (iout,*) 'Too many overlapping confs.',
+ & ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
+ finish=.true.
+ endif
+ else if (.not. finish) then
+C Distribute tasks to processors
+c print *,'Master sending order'
+ call MPI_SEND(12, 1, MPI_INTEGER, is, tag,
+ & CG_COMM, ierr)
+c write (iout,*) '12: tag=',tag
+c print *,'Master sent order to processor',is
+ call MPI_SEND(it, 1, MPI_INTEGER, is, tag,
+ & CG_COMM, ierr)
+c write (iout,*) 'it: tag=',tag
+ call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,
+ & CG_COMM, ierr)
+c write (iout,*) 'eold: tag=',tag
+ call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION,
+ & is, tag,
+ & CG_COMM, ierr)
+c write (iout,*) 'varia: tag=',tag
+ call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION,
+ & is, tag,
+ & CG_COMM, ierr)
+c write (iout,*) 'varold: tag=',tag
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+ imtasks(is)=1
+ imtasks_n=imtasks_n+1
+C End distribution
+ endif ! overlap
+ enddo LOOP3
+
+ flag = .false.
+ LOOP_RECV:do while(.not.flag)
+ do is=1, nodenum-1
+ call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
+ if(flag) then
+ call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,
+ & CG_COMM, status, ierr)
+ call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,
+ & CG_COMM, status, ierr)
+ call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,
+ & CG_COMM, status, ierr)
+ call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,
+ & CG_COMM, status, ierr)
+ call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is,
+ & tag, CG_COMM, status, ierr)
+ call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,
+ & CG_COMM, status, ierr)
+ call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,
+ & CG_COMM, status, ierr)
+ call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,
+ & CG_COMM, status, ierr)
+ i_grnum_d=i_grnum_d+ii_grnum_d
+ i_ennum_d=i_ennum_d+ii_ennum_d
+ neneval = neneval+ii_ennum_d
+ i_hesnum_d=i_hesnum_d+ii_hesnum_d
+ i_minimiz=i_minimiz+1
+ imtasks(is)=0
+ imtasks_n=imtasks_n-1
+ exit
+ endif
+ enddo
+ enddo LOOP_RECV
+
+ if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)')
+ & 'From Worker #',is,' iitt',iitt,
+ & ' Conformation:',ngen,' energy:',etot
+C--------------------------------------------------------------------------
+C... Do Metropolis test
+C--------------------------------------------------------------------------
+ call metropolis(nvar,varia,varold1,etot,eold1,accepted,
+ & similar,EneLower)
+ if(iitt.ne.it.and..not.similar) then
+ call metropolis(nvar,varia,varold,etot,eold,accepted,
+ & similar,EneLower)
+ accepted=enelower
+ endif
+ if(etot.lt.eneglobal)eneglobal=etot
+c if(mod(it,100).eq.0)
+ write(iout,*)'CHUJOJEB ',neneval,eneglobal
+ if (accepted) then
+C Write the accepted conformation.
+ nout=nout+1
+ if (refstr) then
+ call var_to_geom(nvar,varia)
+ call chainbuild
+ call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
+ & nsup,przes,obr,non_conv)
+ rms=dsqrt(rms)
+ call contact(.false.,ncont,icont,co)
+ frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
+ write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
+ & 'RMS deviation from the reference structure:',rms,
+ & ' % of native contacts:',frac*100,' contact order:',co
+ endif ! refstr
+ if (print_mc.gt.0)
+ & write (iout,*) 'Writing new conformation',nout
+ if (print_stat) then
+ call var_to_geom(nvar,varia)
+#if defined(AIX) || defined(PGI)
+ open (istat,file=statname,position='append')
+#else
+ open (istat,file=statname,access='append')
+#endif
+ if (refstr) then
+ write (istat,'(i5,16(1pe14.5))') nout,
+ & (energia(print_order(i)),i=1,nprint_ene),
+ & etot,rms,frac
+ else
+ write (istat,'(i5,16(1pe14.5))') nout,
+ & (energia(print_order(i)),i=1,nprint_ene),etot
+ endif ! refstr
+ close(istat)
+ endif ! print_stat
+C Print internal coordinates.
+ if (print_int) call briefout(nout,etot)
+ nacc=nacc+1
+ nacc_tot=nacc_tot+1
+ if (elowest.gt.etot) elowest=etot
+ moves_acc(MoveType)=moves_acc(MoveType)+1
+ if (MoveType.eq.1) then
+ nbond_acc(nbond)=nbond_acc(nbond)+1
+ endif
+C Check against conformation repetitions.
+ irepet=conf_comp(varia,etot)
+ if (nrepm.gt.maxrepm) then
+ if (print_mc.gt.0)
+ & write (iout,'(a)') 'Too many conformation repetitions.'
+ finish=.true.
+ endif
+C Store the accepted conf. and its energy.
+ eold=etot
+ do i=1,nvar
+ varold(i)=varia(i)
+ enddo
+ if (irepet.eq.0) call zapis(varia,etot)
+C Lower the temperature, if necessary.
+ call cool
+ else
+ ntrial=ntrial+1
+ endif ! accepted
+ 30 continue
+ if(finish.and.imtasks_n.eq.0)exit LOOP2
+ enddo LOOP2 ! accepted
+C Check for time limit.
+ not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
+ if(.not.not_done .or. finish) then
+ if(imtasks_n.gt.0) then
+ not_done=.true.
+ else
+ not_done=.false.
+ endif
+ finish=.true.
+ endif
+ enddo LOOP1 ! not_done
+ runtime=tcpu()
+ if (print_mc.gt.0) then
+ write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
+ call statprint(nacc,nfun,iretcode,etot,elowest)
+ write (iout,'(a)')
+ & 'Statistics of multiple-bond motions. Total motions:'
+ write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
+ write (iout,'(a)') 'Accepted motions:'
+ write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
+ if (it.ge.maxacc)
+ &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
+ endif
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+ do is=1,nodenum-1
+ call MPI_SEND(999, 1, MPI_INTEGER, is, tag,
+ & CG_COMM, ierr)
+ enddo
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine execute_slave(nodeinfo,iprint)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'mpif.h'
+ include 'COMMON.TIME1'
+ include 'COMMON.IOUNITS'
+crc include 'COMMON.DEFORM'
+crc include 'COMMON.DEFORM1'
+crc include 'COMMON.DEFORM2'
+ include 'COMMON.LOCAL'
+ include 'COMMON.VAR'
+ include 'COMMON.INFO'
+ include 'COMMON.MINIM'
+ character*10 nodeinfo
+ double precision x(maxvar),x1(maxvar)
+ nodeinfo='chujwdupe'
+c print *,'Processor:',MyID,' Entering execute_slave'
+ tag=0
+c call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
+c & CG_COMM, ierr)
+
+1001 call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,
+ & CG_COMM, status, ierr)
+c write(iout,*)'12: tag=',tag
+ if(iprint.ge.2)print *, MyID,' recv order ',i_switch
+ if (i_switch.eq.12) then
+ i_grnum_d=0
+ i_ennum_d=0
+ i_hesnum_d=0
+ call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,
+ & CG_COMM, status, ierr)
+c write(iout,*)'12: tag=',tag
+ call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
+ & CG_COMM, status, ierr)
+c write(iout,*)'ener: tag=',tag
+ call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
+ & CG_COMM, status, ierr)
+c write(iout,*)'x: tag=',tag
+ call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
+ & CG_COMM, status, ierr)
+c write(iout,*)'x1: tag=',tag
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+c print *,'calling minimize'
+ call minimize(energyx,x,iretcode,nfun)
+ if(iprint.gt.0)
+ & write(iout,100)'minimized energy = ',energyx,
+ & ' # funeval:',nfun,' iret ',iretcode
+ write(*,100)'minimized energy = ',energyx,
+ & ' # funeval:',nfun,' iret ',iretcode
+ 100 format(a20,f10.5,a12,i5,a6,i2)
+ if(iretcode.eq.10) then
+ do iminrep=2,3
+ if(iprint.gt.1)
+ & write(iout,*)' ... not converged - trying again ',iminrep
+ call minimize(energyx,x,iretcode,nfun)
+ if(iprint.gt.1)
+ & write(iout,*)'minimized energy = ',energyx,
+ & ' # funeval:',nfun,' iret ',iretcode
+ if(iretcode.ne.10)go to 812
+ enddo
+ if(iretcode.eq.10) then
+ if(iprint.gt.1)
+ & write(iout,*)' ... not converged again - giving up'
+ go to 812
+ endif
+ endif
+812 continue
+c print *,'Sending results'
+ call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,
+ & CG_COMM, ierr)
+ call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
+ & CG_COMM, ierr)
+ call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,
+ & CG_COMM, ierr)
+ call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
+ & CG_COMM, ierr)
+ call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
+ & CG_COMM, ierr)
+ call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,
+ & CG_COMM, ierr)
+ call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,
+ & CG_COMM, ierr)
+ call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,
+ & CG_COMM, ierr)
+c print *,'End sending'
+ go to 1001
+ endif
+
+ return
+ end
+#endif
+c------------------------------------------------------------------------------
+ subroutine statprint(it,nfun,iretcode,etot,elowest)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.MCM'
+ if (minim) then
+ write (iout,
+ & '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)')
+ & 'Finished iteration #',it,' energy is',etot,
+ & ' lowest energy:',elowest,
+ & 'SUMSL return code:',iretcode,
+ & ' # of energy evaluations:',neneval,
+ & '# of temperature jumps:',ntherm,
+ & ' # of minima repetitions:',nrepm
+ else
+ write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)')
+ & 'Finished iteration #',it,' energy is',etot,
+ & ' lowest energy:',elowest
+ endif
+ write (iout,'(/4a)')
+ & 'Kind of move ',' total',' accepted',
+ & ' fraction'
+ write (iout,'(58(1h-))')
+ do i=-1,MaxMoveType
+ if (moves(i).eq.0) then
+ fr_mov_i=0.0d0
+ else
+ fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
+ endif
+ write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),
+ & fr_mov_i
+ enddo
+ write (iout,'(a,2i15,f10.5)') 'total ',nmove,nacc_tot,
+ & dfloat(nacc_tot)/dfloat(nmove)
+ write (iout,'(58(1h-))')
+ write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine heat(over)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MCM'
+ include 'COMMON.IOUNITS'
+ logical over
+C Check if there`s a need to increase temperature.
+ if (ntrial.gt.maxtrial) then
+ if (NstepH.gt.0) then
+ if (dabs(Tcur-TMax).lt.1.0D-7) then
+ if (print_mc.gt.0)
+ & write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))')
+ & 'Upper limit of temperature reached. Terminating.'
+ over=.true.
+ Tcur=Tmin
+ else
+ Tcur=Tcur*TstepH
+ if (Tcur.gt.Tmax) Tcur=Tmax
+ betbol=1.0D0/(Rbol*Tcur)
+ if (print_mc.gt.0)
+ & write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
+ & 'System heated up to ',Tcur,' K; BetBol:',betbol
+ ntherm=ntherm+1
+ ntrial=0
+ over=.false.
+ endif
+ else
+ if (print_mc.gt.0)
+ & write (iout,'(a)')
+ & 'Maximum number of trials in a single MCM iteration exceeded.'
+ over=.true.
+ Tcur=Tmin
+ endif
+ else
+ over=.false.
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine cool
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MCM'
+ include 'COMMON.IOUNITS'
+ if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
+ Tcur=Tcur/TstepC
+ if (Tcur.lt.Tmin) Tcur=Tmin
+ betbol=1.0D0/(Rbol*Tcur)
+ if (print_mc.gt.0)
+ & write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
+ & 'System cooled down up to ',Tcur,' K; BetBol:',betbol
+ endif
+ return
+ end
+C---------------------------------------------------------------------------
+ subroutine zapis(varia,etot)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MP
+ include 'mpif.h'
+ include 'COMMON.INFO'
+#endif
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.MCM'
+ include 'COMMON.IOUNITS'
+ integer itemp(maxsave)
+ double precision varia(maxvar)
+ logical lprint
+ lprint=.false.
+ if (lprint) then
+ write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,
+ & ' MaxSave=',MaxSave
+ write (iout,'(a)') 'Current energy and conformation:'
+ write (iout,'(1pe14.5)') etot
+ write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
+ endif
+C Shift the contents of the esave and varsave arrays if filled up.
+ call add2cache(maxvar,maxsave,nsave,nvar,MyID,itemp,
+ & etot,varia,esave,varsave)
+ if (lprint) then
+ write (iout,'(a)') 'Energies and the VarSave array.'
+ do i=1,nsave
+ write (iout,'(i5,1pe14.5)') i,esave(i)
+ write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
+ enddo
+ endif
+ return
+ end
+C---------------------------------------------------------------------------
+ subroutine perturb(error,lprint,MoveType,nbond,max_phi)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ parameter (MMaxSideMove=100)
+ include 'COMMON.MCM'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+crc include 'COMMON.DEFORM1'
+ logical error,lprint,fail
+ integer MoveType,nbond,end_select,ind_side(MMaxSideMove)
+ double precision max_phi
+ double precision psi,gen_psi
+ external iran_num
+ integer iran_num
+ integer ifour
+ data ifour /4/
+ error=.false.
+ lprint=.false.
+C Perturb the conformation according to a randomly selected move.
+ call SelectMove(MoveType)
+c write (iout,*) 'MoveType=',MoveType
+ itrial=0
+ goto (100,200,300,400,500) MoveType
+C------------------------------------------------------------------------------
+C Backbone N-bond move.
+C Select the number of bonds (length of the segment to perturb).
+ 100 continue
+ if (itrial.gt.1000) then
+ write (iout,'(a)') 'Too many attempts at multiple-bond move.'
+ error=.true.
+ return
+ endif
+ bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
+c print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
+c & ' Bond_prob=',Bond_Prob
+ do i=1,nbm-1
+c print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
+ if (bond_prob.ge.sumpro_bond(i) .and.
+ & bond_prob.le.sumpro_bond(i+1)) then
+ nbond=i+1
+ goto 10
+ endif
+ enddo
+ write (iout,'(2a)') 'In PERTURB: Error - number of bonds',
+ & ' to move out of range.'
+ error=.true.
+ return
+ 10 continue
+ if (nwindow.gt.0) then
+C Select the first residue to perturb
+ iwindow=iran_num(1,nwindow)
+ print *,'iwindow=',iwindow
+ iiwin=1
+ do while (winlen(iwindow).lt.nbond)
+ iwindow=iran_num(1,nwindow)
+ iiwin=iiwin+1
+ if (iiwin.gt.1000) then
+ write (iout,'(a)') 'Cannot select moveable residues.'
+ error=.true.
+ return
+ endif
+ enddo
+ nstart=iran_num(winstart(iwindow),winend(iwindow))
+ else
+ nstart = iran_num(koniecl+2,nres-nbond-koniecl)
+cd print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
+cd & ' nstart=',nstart
+ endif
+ psi = gen_psi()
+ if (psi.eq.0.0) then
+ error=.true.
+ return
+ endif
+ if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)')
+ & 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
+cd print *,'nstart=',nstart
+ call bond_move(nbond,nstart,psi,.false.,error)
+ if (error) then
+ write (iout,'(2a)')
+ & 'Could not define reference system in bond_move, ',
+ & 'choosing ahother segment.'
+ itrial=itrial+1
+ goto 100
+ endif
+ nbond_move(nbond)=nbond_move(nbond)+1
+ moves(1)=moves(1)+1
+ nmove=nmove+1
+ return
+C------------------------------------------------------------------------------
+C Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
+C the chain.
+ 200 continue
+ lprint=.true.
+c end_select=iran_num(1,2*koniecl)
+c if (end_select.gt.koniecl) then
+c end_select=nphi-(end_select-koniecl)
+c else
+c end_select=koniecl+3
+c endif
+c if (nwindow.gt.0) then
+c iwin=iran_num(1,nwindow)
+c i1=max0(4,winstart(iwin))
+c i2=min0(winend(imin)+2,nres)
+c end_select=iran_num(i1,i2)
+c else
+c iselect = iran_num(1,nmov_var)
+c jj = 0
+c do i=1,nphi
+c if (isearch_tab(i).eq.1) jj = jj+1
+c if (jj.eq.iselect) then
+c end_select=i+3
+c exit
+c endif
+c enddo
+c endif
+ end_select = iran_num(4,nres)
+ psi=max_phi*gen_psi()
+ if (psi.eq.0.0D0) then
+ error=.true.
+ return
+ endif
+ phi(end_select)=pinorm(phi(end_select)+psi)
+ if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)')
+ & 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',
+ & phi(end_select)*rad2deg
+c if (end_select.gt.3)
+c & theta(end_select-1)=gen_theta(itype(end_select-2),
+c & phi(end_select-1),phi(end_select))
+c if (end_select.lt.nres)
+c & theta(end_select)=gen_theta(itype(end_select-1),
+c & phi(end_select),phi(end_select+1))
+cd print *,'nres=',nres,' end_select=',end_select
+cd print *,'theta',end_select-1,theta(end_select-1)
+cd print *,'theta',end_select,theta(end_select)
+ moves(2)=moves(2)+1
+ nmove=nmove+1
+ lprint=.false.
+ return
+C------------------------------------------------------------------------------
+C Side chain move.
+C Select the number of SCs to perturb.
+ 300 isctry=0
+ 301 nside_move=iran_num(1,MaxSideMove)
+c print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
+C Select the indices.
+ do i=1,nside_move
+ icount=0
+ 111 inds=iran_num(nnt,nct)
+ icount=icount+1
+ if (icount.gt.1000) then
+ write (iout,'(a)')'Error - cannot select side chains to move.'
+ error=.true.
+ return
+ endif
+ if (itype(inds).eq.10) goto 111
+ do j=1,i-1
+ if (inds.eq.ind_side(j)) goto 111
+ enddo
+ do j=1,i-1
+ if (inds.lt.ind_side(j)) then
+ indx=j
+ goto 112
+ endif
+ enddo
+ indx=i
+ 112 do j=i,indx+1,-1
+ ind_side(j)=ind_side(j-1)
+ enddo
+ 113 ind_side(indx)=inds
+ enddo
+C Carry out perturbation.
+ do i=1,nside_move
+ ii=ind_side(i)
+ iti=itype(ii)
+ call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
+ if (fail) then
+ isctry=isctry+1
+ if (isctry.gt.1000) then
+ write (iout,'(a)') 'Too many errors in SC generation.'
+ error=.true.
+ return
+ endif
+ goto 301
+ endif
+ if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)')
+ & 'Side chain ',restyp(iti),ii,' moved to ',
+ & alph(ii)*rad2deg,omeg(ii)*rad2deg
+ enddo
+ moves(3)=moves(3)+1
+ nmove=nmove+1
+ return
+C------------------------------------------------------------------------------
+C THETA move
+ 400 end_select=iran_num(3,nres)
+ theta_new=gen_theta(itype(end_select),phi(end_select),
+ & phi(end_select+1))
+ if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)')
+ & 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,
+ & ' to ',theta_new*rad2deg
+ theta(end_select)=theta_new
+ moves(4)=moves(4)+1
+ nmove=nmove+1
+ return
+C------------------------------------------------------------------------------
+C Error returned from SelectMove.
+ 500 error=.true.
+ return
+ end
+C------------------------------------------------------------------------------
+ subroutine SelectMove(MoveType)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MCM'
+ include 'COMMON.IOUNITS'
+ what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
+ do i=1,MaxMoveType
+ if (what_move.ge.sumpro_type(i-1).and.
+ & what_move.lt.sumpro_type(i)) then
+ MoveType=i
+ return
+ endif
+ enddo
+ write (iout,'(a)')
+ & 'Fatal error in SelectMoveType: cannot select move.'
+ MoveType=MaxMoveType+1
+ return
+ end
+c----------------------------------------------------------------------------
+ double precision function gen_psi()
+ implicit none
+ integer i
+ double precision x,ran_number
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ x=0.0D0
+ do i=1,100
+ x=ran_number(-pi,pi)
+ if (dabs(x).gt.angmin) then
+ gen_psi=x
+ return
+ endif
+ enddo
+ write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
+ gen_psi=0.0D0
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,
+ & enelower)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MCM'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+crc include 'COMMON.DEFORM'
+ double precision ecur,eold,xx,ran_number,bol
+ double precision xcur(n),xold(n)
+ double precision ecut1 ,ecut2 ,tola
+ logical accepted,similar,not_done,enelower
+ logical lprn
+ data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
+! ecut1=-5*enedif
+! ecut2=50*enedif
+! tola=5.0d0
+C Set lprn=.true. for debugging.
+ lprn=.false.
+ if (lprn)
+ &write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
+ similar=.false.
+ enelower=.false.
+ accepted=.false.
+C Check if the conformation is similar.
+ difene=ecur-eold
+ reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
+ if (lprn) then
+ write (iout,*) 'Metropolis'
+ write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
+ endif
+C If energy went down remarkably, we accept the new conformation
+C unconditionally.
+cjp if (reldife.lt.ecut1) then
+ if (difene.lt.ecut1) then
+ accepted=.true.
+ EneLower=.true.
+ if (lprn) write (iout,'(a)')
+ & 'Conformation accepted, because energy has lowered remarkably.'
+! elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
+cjp elseif (reldife.lt.ecut2)
+ elseif (difene.lt.ecut2)
+ & then
+C Reject the conf. if energy has changed insignificantly and there is not
+C much change in conformation.
+ if (lprn)
+ & write (iout,'(2a)') 'Conformation rejected, because it is',
+ & ' similar to the preceding one.'
+ accepted=.false.
+ similar=.true.
+ else
+C Else carry out Metropolis test.
+ EneLower=.false.
+ xx=ran_number(0.0D0,1.0D0)
+ xxh=betbol*difene
+ if (lprn)
+ & write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
+ if (xxh.gt.50.0D0) then
+ bol=0.0D0
+ else
+ bol=exp(-xxh)
+ endif
+ if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
+ if (bol.gt.xx) then
+ accepted=.true.
+ if (lprn) write (iout,'(a)')
+ & 'Conformation accepted, because it passed Metropolis test.'
+ else
+ accepted=.false.
+ if (lprn) write (iout,'(a)')
+ & 'Conformation rejected, because it did not pass Metropolis test.'
+ endif
+ endif
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+ return
+ end
+c------------------------------------------------------------------------------
+ integer function conf_comp(x,ene)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MCM'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ double precision etol , angtol
+ double precision x(maxvar)
+ double precision dif_ang,difa
+ data etol /0.1D0/, angtol /20.0D0/
+ do ii=nsave,1,-1
+c write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
+ if (dabs(ene-esave(ii)).lt.etol) then
+ difa=dif_ang(nphi,x,varsave(1,ii))
+c do i=1,nphi
+c write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
+c & rad2deg*varsave(i,ii)
+c enddo
+c write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
+ if (difa.le.angtol) then
+ if (print_mc.gt.0) then
+ write (iout,'(a,i5,2(a,1pe15.4))')
+ & 'Current conformation matches #',ii,
+ & ' in the store array ene=',ene,' esave=',esave(ii)
+c write (*,'(a,i5,a)') 'Current conformation matches #',ii,
+c & ' in the store array.'
+ endif ! print_mc.gt.0
+ if (print_mc.gt.1) then
+ do i=1,nphi
+ write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
+ & rad2deg*varsave(i,ii)
+ enddo
+ endif ! print_mc.gt.1
+ nrepm=nrepm+1
+ conf_comp=ii
+ return
+ endif
+ endif
+ enddo
+ conf_comp=0
+ return
+ end
+C----------------------------------------------------------------------------
+ double precision function dif_ang(n,x,y)
+ implicit none
+ integer i,n
+ double precision x(n),y(n)
+ double precision w,wa,dif,difa
+ double precision pinorm
+ include 'COMMON.GEO'
+ wa=0.0D0
+ difa=0.0D0
+ do i=1,n
+ dif=dabs(pinorm(y(i)-x(i)))
+ if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
+ w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
+ wa=wa+w
+ difa=difa+dif*dif*w
+ enddo
+ dif_ang=rad2deg*dsqrt(difa/wa)
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,
+ & ecur,xcur,ecache,xcache)
+ implicit none
+ include 'COMMON.GEO'
+ include 'COMMON.IOUNITS'
+ integer n1,n2,ncache,nvar,SourceID,CachSrc(n2)
+ integer i,ii,j
+ double precision ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
+cd write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
+cd write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
+cd write (iout,*) 'Old CACHE array:'
+cd do i=1,ncache
+cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
+cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
+cd enddo
+
+ i=ncache
+ do while (i.gt.0 .and. ecur.lt.ecache(i))
+ i=i-1
+ enddo
+ i=i+1
+cd write (iout,*) 'i=',i,' ncache=',ncache
+ if (ncache.eq.n2) then
+ write (iout,*) 'Cache dimension exceeded',ncache,n2
+ write (iout,*) 'Highest-energy conformation will be removed.'
+ ncache=ncache-1
+ endif
+ do ii=ncache,i,-1
+ ecache(ii+1)=ecache(ii)
+ CachSrc(ii+1)=CachSrc(ii)
+ do j=1,nvar
+ xcache(j,ii+1)=xcache(j,ii)
+ enddo
+ enddo
+ ecache(i)=ecur
+ CachSrc(i)=SourceID
+ do j=1,nvar
+ xcache(j,i)=xcur(j)
+ enddo
+ ncache=ncache+1
+cd write (iout,*) 'New CACHE array:'
+cd do i=1,ncache
+cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
+cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
+cd enddo
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,
+ & xcache)
+ implicit none
+ include 'COMMON.GEO'
+ include 'COMMON.IOUNITS'
+ integer n1,n2,ncache,nvar,CachSrc(n2)
+ integer i,ii,j
+ double precision ecache(n2),xcache(n1,n2)
+
+cd write (iout,*) 'Enter RM_FROM_CACHE'
+cd write (iout,*) 'Old CACHE array:'
+cd do ii=1,ncache
+cd write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
+cd write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
+cd enddo
+
+ do ii=i+1,ncache
+ ecache(ii-1)=ecache(ii)
+ CachSrc(ii-1)=CachSrc(ii)
+ do j=1,nvar
+ xcache(j,ii-1)=xcache(j,ii)
+ enddo
+ enddo
+ ncache=ncache-1
+cd write (iout,*) 'New CACHE array:'
+cd do i=1,ncache
+cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
+cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
+cd enddo
+ return
+ end
projects / unres.git / blobdiff