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,1),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,1), & 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