+++ /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,it)
- 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