X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fmcm.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fmcm.F;h=7f839f4f9c844f6094ad7bf9c71d4dff4ae3af55;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/mcm.F b/source/unres/src_MD-M-newcorr/mcm.F new file mode 100644 index 0000000..7f839f4 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/mcm.F @@ -0,0 +1,1481 @@ + 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) + +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 + kkk=1 + call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk), + & 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