X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fsrc%2Fmcm.F;fp=source%2Funres%2Fsrc_MD%2Fsrc%2Fmcm.F;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=79e567bc82cd1b3beef03539f20227718c1d0d5e;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/src/mcm.F b/source/unres/src_MD/src/mcm.F deleted file mode 100644 index 79e567b..0000000 --- a/source/unres/src_MD/src/mcm.F +++ /dev/null @@ -1,1481 +0,0 @@ - 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