subroutine func1(n,What,x) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" #endif integer n,nf,What double precision x(n) include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.OPTIM" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.PROTFILES" include "COMMON.COMPAR" include "COMMON.PROTNAME" include "COMMON.THERMAL" include "COMMON.TORSION" C Define local variables integer i,ii,iii,ind,jj,j,k,l,ik,jk,iprot,ib,l1,l2,ll1,ll2, & it integer ibatch integer ipass_conf,istart_conf,iend_conf double precision energia(0:max_ene) double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk, & eave_pft(max_ene),emix_pft(max_ene), & elowest_t(2,maxT),esquare_ft, & elowest_aux(2,maxT), & efree_t,emixsq_pft(max_ene), & eneps_mixt(nntyp),eneps_meant(nntyp), & enepsave_ft(nntyp),eneps_mix_ft(nntyp), & eneps_mixsq_ft(nntyp),emean_ft,fac double precision aux double precision tcpu_ini,tcpu logical lprn integer icant external icant character*80 vormat tcpu_ini = tcpu() n_func=n_func+1 c write (iout,*) "FUNC1: Init_Ene ",init_ene call flush(iout) lprn=.true. #ifdef DEBUG write (iout,*) "ELOWEST at the beginning of FUC1" do iprot=1,nprot do ibatch=1,natlike(iprot)+2 do ib=1,nbeta(ibatch,iprot) write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib, & " elowest",elowest(ib,ibatch,iprot) enddo enddo enddo #endif #ifdef MPI call MPI_Bcast( What, 1, MPI_INTEGER, Master, Comm1, IERROR) call MPI_Bcast( cutoffeval, 1, MPI_LOGICAL, Master, Comm1, IERROR) call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION, Master, Comm1, IERROR) if (What.le.0) return #ifdef DEBUG write (iout,*) "Processor",me,me1," What",What write (iout,*) "Processor",me,me1," x",(x(i),i=1,n) call flush(iout) #endif #else #ifdef DEBUG write (iout,*) "x",(x(i),i=1,n) #endif #endif C Convert variables into weights and other UNRES energy parameters. call x2w(n,x) C Calculate the total energies. c write (iout,*) "Processor",me," nprot",nprot c Revise the list of conformations if (What.eq.5) then write (iout,*) "MAKE_LIST called" call flush(iout) call MPI_Bcast( enecut(1), nprot, MPI_DOUBLE_PRECISION, & Master, Comm1, IERROR) call make_list(.true.,.false.,n,x) return endif DO iprot=1,nprot #ifdef DEBUG if (What.eq.2) then open(icsa_rbank,file="corr."//protname(iprot), & access="direct",form="formatted",recl=25+10*natlike(iprot)) vormat='(2i6,2i3,e15.5,' if (natlike(iprot).lt.10) then write(vormat(16:),'(i1,a)') maxdimnat*natlike(iprot),'f7.4/)' else write(vormat(16:),'(i2,a)') maxdimnat*natlike(iprot),'f7.4/)' endif write (iout,*) "format ",vormat call flush(iout) endif #endif #ifdef DEBUG write (iout,*) "Processor",me," iprot",iprot, & " indstart",indstart(me1,iprot)," indend",indend(me1,iprot), & " init_ene",init_ene," mod_fourier",mod_fourier(nloctyp), & " mod_elec",mod_elec," mod_scp",mod_scp," mod_tor",mod_tor, & " mod_sccor",mod_sccor #endif call restore_molinfo(iprot) ii=0 IF (NCHUNK_CONF(IPROT).EQ.1) THEN #ifdef MPI do i=indstart(me1,iprot),indend(me1,iprot) #else do i=1,ntot_work(iprot) #endif c write (iout,*) "i",i," ii",ii," calling ene_eval" call ene_eval(iprot,i,ii) c write (iout,*) "After ene_eval: i",i," ii",ii #ifdef DEBUG if (What.eq.2) then write (icsa_rbank,vormat,rec=i) & i,list_conf(i,iprot), & e_total(i,iprot),((nu(l,k,ii,iprot),l=1,natdim(k,iprot)), & k=1,natlike(iprot)) endif #endif enddo ELSE if (.not.init_ene .and. mod_fourier(nloctyp).eq.0 & .and. .not. mod_elec .and. .not. mod_scp & .and. .not. mod_tor .and. .not. mod_sccor) then open (ientin,file=benefiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec_ene(iprot)) if (natlike(iprot).gt.0) & open (icbase,file=bprotfiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec(iprot)) else open (icbase,file=bprotfiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec(iprot)) c write (iout,*) "lenrec_ene",lenrec_ene(iprot) call flush(iout) open (ientout,file=benefiles(iprot),status="unknown", & form="unformatted",access="direct",recl=lenrec_ene(iprot)) endif ipass_conf=0 #ifdef MPI do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc istart_conf=i iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot)) #else do i=1,ntot_work(iprot),maxstr_proc istart_conf=i iend_conf=min0(i+maxstr_proc-1,ntot(iprot)) #endif ii=0 c c Read the chunk of conformations off a DA scratchfile. c if (.not. init_ene .and. mod_fourier(nloctyp).eq.0 & .and. .not. mod_elec .and. .not. mod_scp & .and. .not. mod_tor .and. .not. mod_sccor) then #ifdef DEBUG write (iout,*) "func1: calling daread_ene",istart_conf, & iend_conf #endif call daread_ene(iprot,istart_conf,iend_conf) if (natlike(iprot).gt.0) then #ifdef DEBUG write (iout,*) "Calling daread_coords",istart_conf, & iend_conf #endif call daread_ccoords(iprot,istart_conf,iend_conf) endif do iii=istart_conf,iend_conf call ene_eval(iprot,iii,ii) enddo #ifdef DEBUG write (iout,*) "func1: calling dawrite_ene",istart_conf, & iend_conf #endif call dawrite_ene(iprot,istart_conf,iend_conf,ientin) else call daread_ccoords(iprot,istart_conf,iend_conf) do iii=istart_conf,iend_conf call ene_eval(iprot,iii,ii) enddo call dawrite_ene(iprot,istart_conf,iend_conf,ientout) endif #ifdef DEBUG if (What.eq.2) then do iii=istart_conf,iend_conf write (icsa_rbank,vormat,rec=iii) & iii,list_conf(iii,iprot), & e_total(iii,iprot), & ((nu(l,k,iii-istart_conf+1,iprot),l=1,natdim(k,iprot)), & k=1,natlike(iprot)) enddo endif #endif enddo if (.not.init_ene .and. mod_fourier(nloctyp).eq.0 & .and. .not. mod_elec .and. .not. mod_scp & .and. .not. mod_tor .and. .not. mod_sccor) then close (ientin) if (natlike(iprot).gt.0) close(icbase) else close (icbase) endif if (init_ene .or. mod_fourier(nloctyp).gt.0 .or. mod_elec & .or. mod_scp .or. mod_tor .or. mod_sccor) close(ientout) ENDIF c print *,"Processor",me,me1," finished computing energies" c call flush(iout) #ifdef MPI if (What.eq.2 .or. What.eq.3) then call MPI_Gatherv(e_total(indstart(me1,iprot),iprot), & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) call MPI_Gatherv(entfac(indstart(me1,iprot),iprot), & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) call MPI_Gatherv(eini(indstart(me1,iprot),iprot), & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot), & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) endif if (What.eq.3) then call MPI_Allgatherv(e_total(indstart(me1,iprot),iprot), & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION, & Comm1, IERROR) endif #endif c Mean and lowest energies of structural classes IF (NCHUNK_CONF(IPROT).EQ.1) THEN #ifdef MPI do i=1,ntot_work(iprot) i2ii(i,iprot)=0 enddo ii=0 do i=indstart(me1,iprot),indend(me1,iprot) ii=ii+1 i2ii(i,iprot)=ii enddo istart_conf=indstart(me1,iprot) iend_conf=indend(me1,iprot) #else do i=1,ntot_work(iprot) i2ii(i,iprot)=i endif istart_conf=1 iend_conf=ntot_work(iprot) #endif #ifdef DEBUG do i=1,ntot_work(iprot) write (iout,*) "i",i," i2ii",i2ii(i,iprot) enddo call flush(iout) #endif open (ientin,file=benefiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec_ene(iprot)) c call daread_ene(iprot,istart_conf,iend_conf) call emin_search(iprot) ELSE open (ientin,file=benefiles(iprot),status="old", & form="unformatted",access="direct",recl=lenrec_ene(iprot)) ipass_conf=0 #ifdef MPI do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot)) #else do istart_conf=1,ntot_work(iprot),maxstr_proc iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot)) #endif c c Read the chunk of energies and derivatives off a DA scratchfile. c ipass_conf=ipass_conf+1 do i=1,ntot_work(iprot) i2ii(i,iprot)=0 enddo ii=0 do i=istart_conf,iend_conf ii=ii+1 i2ii(i,iprot)=ii enddo #ifdef DEBUG write (iout,*) "ipass_conf",ipass_conf, & " istart_conf",istart_conf," iend_conf",iend_conf do i=1,ntot_work(iprot) write (iout,*) "i",i," i2ii",i2ii(i,iprot) enddo call flush(iout) #endif call daread_ene(iprot,istart_conf,iend_conf) call emin_search(iprot) enddo close(ientin) ENDIF ENDDO ! iprot #ifdef MPI c Complete the calculation of the lowest energies over all classes and c distribute the values to all procs do iprot=1,nprot #ifdef DEBUG write (iout,*) "Processor",me,me1," calling MPI_Reduce: 2" do k=1,natlike(iprot)+2 write (iout,*) "iprot",iprot," k",k do ib=1,nbeta(k,iprot) write (iout,*) "ib",ib write (iout,'(7hELOWEST,10e15.3)') & elowest(ib,k,iprot) enddo enddo call flush(iout) #endif do ibatch=1,natlike(iprot)+2 do ib=1,nbeta(ibatch,iprot) elowest_aux(1,ib)=elowest(ib,ibatch,iprot) elowest_aux(2,ib)=ind_lowest(ib,ibatch,iprot) enddo call MPI_Allreduce(elowest_aux(1,1),elowest_t(1,1), & nbeta(ibatch,iprot), & MPI_2DOUBLE_PRECISION, MPI_MINLOC, Comm1, IERROR) #ifdef DEBUG do ib=1,nbeta(ibatch,iprot) write (iout,*) "beta=",betaT(ib,ibatch,iprot) write (iout,'(9helowest_t,10f15.3)') & elowest_t(1,ib),elowest_t(2,ib) enddo write (iout,*) "Processor",me,me1," finished MPI_Reduce: 2" #endif do ib=1,nbeta(ibatch,iprot) c write (iout,*) "IB",ib," ibatch",ibatch," iprot",iprot call flush(iout) elowest(ib,ibatch,iprot)=elowest_t(1,ib) c write (iout,*) "elowest",ib,ibatch,iprot, c & elowest(ib,ibatch,iprot) call flush(iout) ind_lowest(ib,ibatch,iprot)=elowest_t(2,ib) enddo enddo ! ibatc enddo ! iprot #endif do iprot=1,nprot call averages(iprot) c write (iout,*) "After averages, calling Thermal" c call flush(iout) if (What.eq.2) call thermal(iprot) c write (iout,*) "After thermal" c call flush(iout) enddo ! iprot #ifdef MPI if (me1.eq.Master) then #endif do iprot=1,nprot #ifdef DEBUG write (iout,*) "iprot",iprot," elowest", & ((elowest(ib,k,iprot),ib=1,nbeta(k,iprot)),k=1,natlike(iprot)+2) call flush(iout) #endif do ib=1,nbeta(2,iprot) fac=betaT(ib,2,iprot) c 6/15/05 AL Heat capacity zvtot(ib,iprot)=(esquare_ftot(ib,iprot)- & emean_ftot(ib,iprot)**2)*fac*fac*Rgas & -ebis_ftot(ib,iprot)+heatbase(iprot) #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib, & " emeantot",emean_ftot(ib,iprot), & " esquaretot",esquare_ftot(ib,iprot), & " ebistot",ebis_ftot(ib,iprot), & " zvtot",zvtot(ib,iprot) #endif #ifdef DEBUG write (iout,*) "beta", & betaT(ib,2,iprot), & " efree",efree(ib,2,iprot), & " emean",emean_ftot(ib,iprot), & " variance",esquare_ftot(ib,iprot) call flush(iout) #endif c 4/13/04 AL Correlation coefficients between energy and nativelikeness c in the respective classes. #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib, & " ii",ii write (iout,*) "emean",emean_ftot(ib,iprot), & " esquare",esquare_ftot(ib,iprot) #endif enddo enddo #ifdef DEBUG write (iout,*) & "Averages of nativelikeness and energy-nativelikeness", & " correlation coefficients:" do iprot=1,nprot do i=1,natlike(iprot) do ib=1,nbeta(i+1,iprot) write (iout,*) "protein",iprot," property",i, & " temperature",ib," natdim",natdim(i,iprot)," value(s)", & (nuave(j,ib,i,iprot),j=1,natdim(i,iprot)) enddo enddo enddo call flush(iout) #endif lprn=.false. #ifdef MPI endif #endif #ifdef DEBUG write (iout,*) "Calling CUTOFF_VIOLATION" call flush(iout) #endif c write (iout,*) "CUTOFFEVAL",cutoffeval cutoffviol=.false. if (cutoffeval) call cutoff_violation t_func = t_func + tcpu() - tcpu_ini #ifdef DEBUG write (iout,*) "Exitting FUNC1" call flush(iout) #endif return end c------------------------------------------------------------------------------- subroutine targetfun(n,x,nf,f,uiparm,ufparm) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,Errcode,status(MPI_STATUS_SIZE) include "COMMON.MPI" #endif integer n,nf double precision x(n) include "COMMON.WEIGHTS" include "COMMON.IOUNITS" include "COMMON.DERIV" include "COMMON.VMCPAR" include "COMMON.CLASSES" include "COMMON.GEO" include "COMMON.ENERGIES" include "COMMON.WEIGHTDER" include "COMMON.OPTIM" include "COMMON.TIME1" include "COMMON.COMPAR" double precision ff common /sru/ ff integer i,ii,it,inat,j,k,l,kk,iprot,ib logical bs external ufparm double precision urparm integer uiparm double precision f,kara,viol,gnmr,gnmr1,aux logical lprn lprn = .true. #ifdef DEBUG write (iout,*) me,me1,"Calling TARGETFUN" #endif call write_restart(n,x) nfl=nf call x2w(n,x) f=0.0d0 call func1(n,1,x) f = 0.0d0 c Maximum likelihood contribution do iprot=1,nprot do iT=1,nbeta(1,iprot) f = f + weilik(iT,iprot)*sumlik(iT,iprot) enddo enddo #ifdef DEBUG write (iout,*)"target_fun: sumlik" do iprot=nprot write (iout,*) (sumlik(ii,iprot),ii=1,nbeta_l(iprot)) enddo write (iout,*) "f before Cv =",f call flush(iout) #endif do iprot=1,nprot c 6/15/05 AL Contribution from Cv do ib=1,nbeta(2,iprot) f=f+weiCv(ib,iprot)*(zvtot(ib,iprot)-target_cv(ib,iprot))**2 enddo c write (iout,*) "target_fun from CV: f=",f enddo c write (iout,*) "target_fun: f=",f c 8/2/2013 Contribution from conformation-dependent averages do iprot=1,nprot do inat=1,natlike(iprot) do ib=1,nbeta(inat+2,iprot) do k=1,natdim(inat,iprot) f = f+weinat(k,ib,inat,iprot)* & (nuave(k,ib,inat,iprot)-nuexp(k,ib,inat,iprot))**2 #ifdef DEBUG write (iout,*) "iprot",iprot," inat",inat, & " ib",ib," k=",k," nuave",nuave(k,ib,inat,iproti), & " nuexp",nuexp(k,ib,inat,iprot), & " weight",weinat(k,ib,inat,iprot) #endif enddo ! k enddo ! ib enddo ! inat enddo ! iprot c write (iout,*) "target_fun after adding nativelikeness: f=",f if (mode.lt.4) then c add constraints on weights kara=0.0d0 do i=1,n kara=kara+gnmr1(x(i),x_low(i),x_up(i)) enddo f=f+1.0d16*kara endif ff = f c write (iout,*) "target_fun: f=",f,1.0d16*kara call flush(iout) return end c----------------------------------------------------------------------------- subroutine targetgrad(n,x,nf,g,uiparm,rdum,ufparm) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "COMMON.TIME1" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,status(MPI_STATUS_SIZE) #endif integer n,nf double precision x(n), & g(n),gg(max_paropt) include "COMMON.WEIGHTS" include "COMMON.ENERGIES" include "COMMON.CLASSES" include "COMMON.IOUNITS" include "COMMON.DERIV" include "COMMON.COMPAR" include "COMMON.WEIGHTDER" include "COMMON.OPTIM" include "COMMON.VMCPAR" include "COMMON.GEO" include "COMMON.NAMES" include "COMMON.VAR" double precision ff common /sru/ ff #ifdef MPI include "COMMON.MPI" #endif double precision xi,delta /1.0d-6/ integer uiparm double precision urparm external ufparm double precision zewn(MaxT,maxprot),aux,zscder double precision zewnvtot(maxT,maxprot) double precision sumlik0(maxT,maxprot),zv0tot(maxT,maxprot), & nuave0(maxdimnat,maxT,maxnatlike,maxprot) double precision geps(nntyp),gtor(maxtor_var) double precision rdum double precision fac integer ll1,ll2,nl1,nl2,nn integer i,ii,ind,inat,j,k,kk,iprot,ib,n_weight integer iti,itj,jj,icant integer l,ll logical lprn double precision gnmrprim,gnmr1prim double precision tcpu_ini,tcpu logical in_sc_pair_list external in_sc_pair_list lprn=.false. tcpu_ini = tcpu() n_grad = n_grad+1 #ifdef MPI #ifdef DEBUG write (iout,*) "Processor",me,me1," calling TARGETGRAD" call flush(iout) #endif #else #ifdef DEBUG write (iout,*) "Calling TARGETGRAD" call flush(iout) #endif #endif #ifdef DEBUG write (iout,*) "x",(x(i),i=1,n) call flush(iout) #endif icg=mod(nf,2)+1 call x2w(n,x) call func1(n,1,x) do iprot=1,nprot do ib=1,nbeta(2,iprot) c 6/15/05 AL Contributions from Cv zewnvtot(ib,iprot)=2*weiCv(ib,iprot)* & (zvtot(ib,iprot)-target_cv(ib,iprot)) enddo enddo do k=1,n g(k)=0.0d0 enddo if (mod_side) then do k=1,nntyp geps(k)=0.0d0 enddo endif if (mod_tor .or. mod_sccor) then do k=1,maxtor_var gtor(k)=0.0d0 enddo endif c Maximum likelihood gradient do iprot=1,nprot do ib=1,nbeta(1,iprot) kk=0 do k=1,n_ene if (imask(k).gt.0) then kk=kk+1 g(kk)=g(kk)+weilik(ib,iprot)*sumlikder(kk,ib,iprot) endif enddo do k=1,nntyp geps(k)=geps(k)+weilik(ib,iprot)*sumlikeps(k,ib,iprot) enddo do k=1,ntor_var gtor(k)=gtor(k)+weilik(ib,iprot)*sumliktor(k,ib,iprot) enddo enddo enddo #ifdef DEBUG write (iout,*) "gtor" do k=1,ntor_var write(iout,*) k,gtor(k) enddo #endif call zderiv c Heat-capacity gradient do iprot=1,nprot do ib=1,nbeta(2,iprot) kk=0 do k=1,n_ene if (imask(k).gt.0) then kk=kk+1 g(kk)=g(kk)+zvdertot(kk,ib,iprot)*zewnvtot(ib,iprot) endif enddo do k=1,nntyp geps(k)=geps(k)+zewnvtot(ib,iprot)*zvepstot(k,ib,iprot) enddo do k=1,ntor_var gtor(k)=gtor(k)+zewnvtot(ib,iprot)*zvtortot(k,ib,iprot) enddo enddo enddo #ifdef DEBUG write (iout,*) "gtor" do k=1,ntor_var write(iout,*) k,gtor(k) enddo write (iout,*) "grad: g",(g(i),i=1,n) call flush(iout) #endif c Derivatives of conformational-dependent properties do iprot=1,nprot do inat=1,natlike(iprot) do ib=1,nbeta(inat+2,iprot) call propderiv(ib,inat,iprot) do l=1,natdim(inat,iprot) aux=2*weinat(l,ib,inat,iprot)* & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot)) #ifdef DEBUG write (iout,*) "iprot",iprot," ib",ib, & " inat",inat," l",l, & " nuave",nuave(l,ib,inat,iprot), & " weight",weinat(ib,inat,iprot), & " zewnR",aux call flush(iout) #endif kk=0 do k=1,n_ene if (imask(k).gt.0) then kk=kk+1 g(kk)=g(kk)+aux*rder(kk,l) endif enddo ! k do k=1,nntyp geps(k)=geps(k)+aux*reps(k,l) enddo ! k do k=1,ntor_var gtor(k)=gtor(k)+aux*rtor(k,l) enddo ! k enddo ! l enddo ! ib enddo enddo ! iprot #ifdef DEBUG write (iout,*) "gtor" do k=1,ntor_var write(iout,*) k,gtor(k) enddo #endif c Compute numerically the remaining part of the gradient n_weight=0 do i=1,n_ene if (imask(i).gt.0) n_weight=n_weight+1 enddo ii=n_weight c Sum up the gradient in SC epsilons, if needed if (mod_side) then do i=1,nsingle_sc ii=ii+1 iti=ityp_ssc(i) g(ii)=g(ii)+geps(icant(iti,iti)) do j=1,iti-1 if (.not. in_sc_pair_list(iti,j)) then g(ii)=g(ii)+0.5d0*geps(icant(iti,j)) endif enddo enddo do i=1,npair_sc ii=ii+1 iti=ityp_psc(1,i) itj=ityp_psc(2,i) g(ii)=geps(icant(iti,itj)) enddo endif c write (iout,*) "n_ene",n_ene if (mod_tor .or. mod_sccor) then do i=1,ntor_var ii=ii+1 g(ii)=gtor(i) enddo endif n_weight=ii #ifdef DEBUG write (iout,*) "n",n," n_weight",n_weight write (iout,*) "After SC: grad:" do i=1,n write (iout,*) i,g(i) enddo call flush(iout) #endif call func1(n,1,x) do iprot=1,nprot do ib=1,nbeta(1,iprot) sumlik0(ib,iprot)=sumlik(ib,iprot) enddo enddo do iprot=1,nprot do ib=1,nbeta(2,iprot) zv0tot(ib,iprot)=zvtot(ib,iprot) enddo enddo do iprot=1,nprot do i=1,natlike(iprot) do ib=1,nbeta(i+2,iprot) do k=1,natdim(i,iprot) nuave0(k,ib,i,iprot)=nuave(k,ib,i,iprot) enddo enddo enddo enddo nn=n do iprot=1,nprot if (nbeta(2,iprot).gt.0) nn=nn-1 enddo ii=nn do iprot=1,nprot if (nbeta(2,iprot).gt.0) then ii=ii+1 do ib=1,nbeta(2,iprot) g(ii)=g(ii)+zewnvtot(ib,iprot) enddo endif enddo #ifdef DEBUG write (iout,*) "n",n," n_weight",n_weight write (iout,*) "After heat grad:" do i=1,n write (iout,*) i,g(i) enddo call flush(iout) #endif c write (iout,*) "n_weight",n_weight," n",n do i=n_weight+1,nn xi=x(i) x(i)=xi+delta c write (iout,*) "i",i," x",x(i),xi call x2w(n,x) call func1(n,1,x) c Maximum likelihood do iprot=1,nprot do j=1,nbeta(1,iprot) g(i)=g(i)+weilik(ib,iprot)*(sumlik(j,iprot) & -sumlik0(j,iprot))/delta enddo enddo #ifdef DEBUG write (iout,*) "init numerical derivatives" call flush(iout) #endif c Contribution from Cv do iprot=1,nprot do ib=1,nbeta(2,iprot) g(i)=g(i)+zewnvtot(ib,iprot)* & (zvtot(ib,iprot)-zv0tot(ib,iprot))/delta enddo enddo c Contribution from average nativelikenss and correlation coefficients c between energy and nativelikeness do iprot=1,nprot do ii=1,natlike(iprot) do ib=1,nbeta(ii+2,iprot) do k=1,natdim(i,iprot) aux=2*weinat(k,ib,inat,iprot)* & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot)) g(i)=g(i)+aux*(nuave(k,ib,ii,iprot)- & nuave0(k,ib,ii,iprot))/delta enddo ! k enddo enddo enddo x(i)=xi enddo #ifdef DEBUG write (iout,*) "n",n," n_weight",n_weight write (iout,*) "After num grad:" do k=1,n write (iout,*) k,g(k) enddo call flush(iout) #endif call func1(n,1,x) #ifdef DEBUG write (iout,*) "finish numerical derivatives" call flush(iout) #endif c c 1/10/2007 AL Correction: ENDIF...IF insterted to prevent unmatched c SEND is I_LOCAL_CHECK.GT.0 c #ifdef DEBUG write (iout,*) "mode",mode #endif if (mode.gt.3) then c transform the gradient to the x space do i=1,n g(i)=g(i)*(x_up(i)-x_low(i))/(pi*(1.0d0+(x(i) & -0.5d0*(x_up(i)+x_low(i)))**2)) enddo else c add constraints on weights c write (iout,*) "Gradient in constraints" do i=1,n g(i)=g(i)+1.0d16*gnmr1prim(x(i),x_low(i),x_up(i)) #ifdef DEBUG write (iout,*) i,x(i),x_low(i),x_up(i), & gnmr1prim(x(i),x_low(i),x_up(i)) #endif enddo endif #ifdef DEBUG write (iout,*) "grad" do i=1,n write (iout,*) i,g(i) enddo #endif t_grad = t_grad + tcpu() - tcpu_ini #ifdef DEBUG write (iout,*) "Exit targetgrad" call flush(iout) #endif return end c------------------------------------------------------------------------------ subroutine matout(ndim,n,a) implicit real*8 (a-h,o-z) include "DIMENSIONS.ZSCOPT" include "COMMON.IOUNITS" double precision a(ndim) character*6 line6 / '------' / character*12 line12 / '------------' / dimension b(8) do 1 i=1,n,6 nlim=min0(i+5,n) write (iout,1000) (k,k=i,nlim) write (iout,1020) line6,(line12,k=i,nlim) 1000 format (/7x,6(i7,5x)) 1020 format (a6,6a12) do 2 j=i,n jlim=min0(j,nlim) jlim1=jlim-i+1 do 3 k=i,jlim 3 b(k-i+1)=a(icant(i,j)) write (iout,1010) j,(b(k),k=1,jlim1) 2 continue 1 continue 1010 format (i3,3x,6(1pe11.4,1x)) return end c------------------------------------------------------------------------------ subroutine zderiv implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,status(MPI_STATUS_SIZE) #endif integer n,nf include "COMMON.WEIGHTS" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.DERIV" include "COMMON.WEIGHTDER" include "COMMON.VMCPAR" include "COMMON.GEO" include "COMMON.NAMES" include "COMMON.VAR" include "COMMON.CLASSES" include "COMMON.THERMAL" #ifdef MPI include "COMMON.MPI" #endif include "COMMON.TIME1" double precision aux,zscder double precision fac,fac2 double precision efj(2),emj(2),varj(2),edecum,edecum1,ebisdecum, & e2decum double precision eavpfj(max_ene,2),emixpfj(max_ene,2), & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2), & enepsmixsqfj(nntyp,2),eavpfjprim(max_ene,2), & entoravefj(0:3,maxtor_var,2),entoravefjprim(0:3,maxtor_var,2), & entormixfj(0:3,maxtor_var,2),entormixsqfj(0:3,maxtor_var,2) integer i,ii,j,k,kk,iprot,ib integer iti,itj,jj,icant integer l,ll logical lprn lprn=.false. do iprot=1,nprot C Heat capacity do ib=1,nbeta(2,iprot) fac = betaT(ib,2,iprot) fac2 = fac*fac kk=0 do k=1,n_ene if (imask(k).gt.0) then kk=kk+1 ebisdecum=emix_pfbistot(k,ib,iprot)- & eave_pftot(k,ib,iprot) & *ebis_ftot(ib,iprot) edecum=emix_pfprimtot(k,ib,iprot) & -eave_pfprimtot(k,ib,iprot) & *emean_ftot(ib,iprot) edecum1=emix_pftot(k,ib,iprot) & -eave_pftot(k,ib,iprot) & *emean_ftot(ib,iprot) e2decum=emixsq_pftot(k,ib,iprot) & -eave_pftot(k,ib,iprot) & *esquare_ftot(ib,iprot) zvdertot(kk,ib,iprot)=Rgas*fac2*(2*edecum & -fac*(e2decum-2*emean_ftot(ib,iprot) & *edecum1)) & -eave_pfbistot(k,ib,iprot) & +fac*ebisdecum #ifdef DEBUG write (iout,*) "k=",k, & "eave_pftot",eave_pftot(k,ib,iprot), & " emix_pftot",emix_pftot(k,ib,iprot), & " emix_pfprimtot",emix_pfprimtot(k,ib,iprot), & " emix_pfbistot",emix_pfbistot(k,ib,iprot), & " eave_pfprimtot",eave_pfprimtot(k,ib,iprot), & " eave_pfbistot",eave_pfbistot(k,ib,iprot), & " emixsq_pftot",emixsq_pftot(k,ib,iprot), & " ebis_ftot",ebis_ftot(ib,iprot), & " emean_ftot",emean_ftot(ib,iprot) write (iout,*) "k=",k," edecum",edecum, & " edecum1",edecum1, & " e2decum",e2decum," ebisdecum",ebisdecum, & " zvdertot", & zvdertot(kk,ib,iprot) #endif endif enddo if (mod_side) then do k=1,nntyp ebisdecum=eneps_mix_fbistot(k,ib,iprot)- & enepsave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot) edecum=eneps_mix_ftot(k,ib,iprot) & -enepsave_ftot(k,ib,iprot)*emean_ftot(ib,iprot) e2decum=eneps_mixsq_ftot(k,ib,iprot) & -enepsave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot) zvepstot(k,ib,iprot)=ww(1)*(Rgas*fac2*(2*edecum & -fac*(e2decum-2*emean_ftot(ib,iprot) & *edecum))+fac*ebisdecum) #ifdef DEBUG write (iout,*) ib,k,"ebisdecum",ebisdecum, & " edecum",edecum," e2decum",e2decum," zvepstot", & zvepstot(k,ibb,iprot) write (iout,*) k, & " enepsave_ftot",enepsave_ftot(k,ib,iprot), & " eneps_mix_ftot",eneps_mix_ftot(k,ib,iprot), & " eneps_mix_fbistot",eneps_mix_fbistot(k,ib,iprot), & " eneps_mixsq_ftot",eneps_mixsq_ftot(k,ib,iprot), & " emean_ftot",emean_ftot(ib,iprot), & " ebis_ftot",ebis_ftot(ib,iprot) #endif enddo endif if (mod_tor .or. mod_sccor) then do k=1,ntor_var ebisdecum=entor_mix_fbistot(k,ib,iprot)- & entorave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot) edecum=entor_mix_fprimtot(k,ib,iprot) & -entorave_fprimtot(k,ib,iprot)*emean_ftot(ib,iprot) edecum1=entor_mix_ftot(k,ib,iprot) & -entorave_ftot(k,ib,iprot)*emean_ftot(ib,iprot) e2decum=entor_mixsq_ftot(k,ib,iprot) & -entorave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot) zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum & -fac*(e2decum-2*emean_ftot(ib,iprot) & *edecum1))-entorave_fbistot(k,ib,iprot) & +fac*ebisdecum #ifdef DEBUG write (iout,*) k, & " entorave_ftot",entorave_ftot(k,ib,iprot), & " entor_mix_ftot",entor_mix_ftot(k,ib,iprot), & " entor_mix_fprimtot",entor_mix_fprimtot(k,ib,iprot), & " entor_mix_fbistot",entor_mix_fbistot(k,ib,iprot), & " entorave_fprimtot",entorave_fprimtot(k,ib,iprot), & " entorave_fbistot",entorave_fbistot(k,ib,iprot), & " entor_mixsq_ftot",entor_mixsq_ftot(k,ib,iprot), & " emean_ftot",emean_ftot(ib,iprot), & " ebis_ftot",ebis_ftot(ib,iprot) write (iout,*) ib,k," ebisdecum",ebisdecum, & " edecum",edecum," edecum1",edecum1," e2decum",e2decum, & " zvtortot",zvtortot(k,ibb,iprot) #endif enddo endif #ifdef DEBUG write (iout,*) "iprot",iprot, & " ib",ib write (iout,*) "zvder", & (zvdertot(k,ibb,iprot),k=1,kk) #endif enddo ! ib enddo ! iprot return end c------------------------------------------------------------------------------ subroutine propderiv(ib,inat,iprot) c Compute the derivatives of conformation-dependent averages in force-field c parameters. implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,status(MPI_STATUS_SIZE) #endif integer n,nf,inat include "COMMON.WEIGHTS" include "COMMON.ENERGIES" include "COMMON.CLASSES" include "COMMON.IOUNITS" include "COMMON.DERIV" include "COMMON.WEIGHTDER" include "COMMON.VMCPAR" include "COMMON.GEO" include "COMMON.OPTIM" include "COMMON.NAMES" include "COMMON.COMPAR" #ifdef MPI include "COMMON.MPI" #endif include "COMMON.TIME1" include "COMMON.VAR" double precision aux,zscder double precision fac double precision efj(2),emj(2),varj(2) double precision eavpfj(max_ene,2),emixpfj(max_ene,2), & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2), & enepsmixsqfj(nntyp,2) integer ll1,nl1 integer i,ii,ind,j,k,kk,iprot,ib,ibb,l1 integer iti,itj,jj,icant integer l,ll logical lprn c write (iout,*) "Entered PROPDERIV" lprn=.false. fac = betaT(ib,inat+2,iprot) c write (iout,*) "derivatives: iprot",iprot, c & " ib",ib," l",l do l=1,natdim(inat,iprot) kk=0 do k=1,n_ene if (imask(k).gt.0) then kk=kk+1 rder(kk,l)=-fac*(nu_pf(k,l,ib,inat,iprot) & -nuave(l,ib,inat,iprot)* & eave_nat_pftot(k,ib,inat,iprot)) endif enddo enddo if (mod_side) then do l=1,natdim(inat,iprot) do k=1,nntyp reps(k,l)=-ww(1)*fac*(nuepsave_f(k,l,ib,inat,iprot)- & nuave(l,ib,inat,iprot)*enepsave_nat_ftot(k,ib,inat,iprot)) enddo enddo endif if (mod_tor .or. mod_sccor) then do l=1,natdim(inat,iprot) do k=1,ntor_var rtor(k,l)=-fac*(nutorave_f(k,l,ib,inat,iprot)- & nuave(l,inat,ib,iprot)*entorave_nat_ftot(k,ib,inat,iprot)) enddo enddo endif return end c------------------------------------------------------------------------------- subroutine ene_eval(iprot,i,ii) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" #ifdef MPI include "mpif.h" integer IERROR,ErrCode include "COMMON.MPI" #endif include "COMMON.CONTROL" include "COMMON.COMPAR" include "COMMON.WEIGHTS" include "COMMON.WEIGHTDER" include "COMMON.CLASSES" include "COMMON.ENERGIES" include "COMMON.IOUNITS" include "COMMON.VMCPAR" include "COMMON.NAMES" include "COMMON.INTERACT" include "COMMON.TIME1" include "COMMON.CHAIN" include "COMMON.VAR" include "COMMON.GEO" include "COMMON.SCCOR" include "COMMON.TORSION" C Define local variables integer i,ii,ij,jj,j,k,ik,jk,l,ll,lll,iprot,iblock double precision energia(0:max_ene) double precision etoti,enepsjk double precision ftune_eps external ftune_eps logical lprn integer icant external icant lprn=.true. ii = ii+1 #ifdef DEBUG write (iout,*) "ene_eval: iprot",iprot," i",i," ii",ii #endif if (init_ene .or. (mod_fourier(nloctyp).gt.0 & .or. mod_elec .or. mod_scp .or. mod_tor .or. mod_sccor) ) then #ifdef DEBUG write (iout,*) "FUNC1: doing UNRES" write (iout,*) "Protein",iprot," i",i," ii",ii," nres",nres call flush(iout) #endif call restore_ccoords(iprot,ii) call int_from_cart1(.false.) #ifdef DEBUG write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct) write (iout,'(8f10.4)') (vbld(k),k=2,nres) write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) c call intout c call cartprint #endif call etotal(energia(0)) c if (natlike(iprot).gt.0) c & call natsimil(i,ii,iprot,.false.) e_total(i,iprot)=energia(0) do j=1,n_ene enetb(ii,j,iprot)=energia(j) enddo if (mod_side) then write (iout,*) "FUN nntyp",nntyp," eneps",eneps_temp(1,45) do j=1,nntyp do k=1,2 eneps(k,j,ii,iprot)=eneps_temp(k,j) enddo enddo #ifdef DEBUG write (iout,'(2i5,18(1pe12.4))') i,ii, & (enetb(ii,j,iprot),j=1,n_ene) write(iout,'(4hENER,i10,3f15.3,2i10,i5)') & i,energia(0),eini(i,iprot),entfac(i,iprot) #endif c#ifdef DEBUG c write (iout,*) "Matrix eneps" c do j=1,ntyp c do k=1,j c write (iout,*) j,k,icant(j,k), c & eneps_temp(1,icant(j,k)), c & eneps_temp(2,icant(j,k)) c enddo c enddo c call flush(iout) c#endif endif if (mod_tor .or. mod_sccor) then c write (iout,*) "etor",enetb(ii,13,iprot),enetb(ii,19,iprot) k=0 do iblock=1,2 do ij=-ntortyp+1,ntortyp-1 do jj=-ntortyp+1,ntortyp-1 if (mask_tor(0,jj,ij,iblock).eq.1) then k=k+1 entor(k,ii,iprot)=etor_temp(0,0,jj,ij,iblock) c write (iout,*) "etor_temp",restyp(itortyp(j)), c & restyp(itortyp(i)), c & k,etor_temp(ji,ii) else if (mask_tor(0,jj,ij,iblock).gt.1) then if (tor_mode.eq.0) then do l=1,nterm(jj,ij,iblock) k=k+1 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock) enddo if (mask_tor(0,jj,ij,iblock).gt.2) then do l=nterm(jj,ij,iblock)+1,2*nterm(jj,ij,iblock) k=k+1 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock) enddo endif else if (tor_mode.eq.1) then do l=1,nterm_kcc(jj,ij) do ll=1,nterm_kcc_Tb(jj,ij) do lll=1,nterm_kcc_Tb(jj,ij) k=k+1 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj) enddo ! lll enddo ! ll enddo if (mask_tor(0,jj,ij,1).gt.2) then do l=nterm_kcc(jj,ij)+1,2*nterm_kcc(jj,ij) do ll=1,nterm_kcc_Tb(jj,ij) do lll=1,nterm_kcc_Tb(jj,ij) k=k+1 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj) enddo enddo enddo endif endif endif enddo enddo enddo do ij=-nsccortyp+1,nsccortyp-1 do jj=-nsccortyp+1,nsccortyp-1 do ll=1,3 if (mask_tor(ll,jj,ij,1).eq.1) then k=k+1 entor(k,ii,iprot)=etor_temp(0,ll,jj,ij,1) c write (iout,*) "etor_temp",restyp(isccortyp(jj)), c & restyp(isccortyp(ij)), c & k,etor_temp(jj,ij) else if (mask_tor(ll,jj,ij,1).gt.1) then do l=1,nterm_sccor(jj,ij) k=k+1 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1) enddo do l=nterm_sccor(jj,ij)+1,2*nterm_sccor(jj,ij) k=k+1 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1) enddo endif enddo enddo enddo endif #ifdef DEBUG c write (iout,*) "Conformation",i c write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres) c write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct) c write (iout,'(8f10.4)') (vbld(k),k=2,nres) c write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct) c write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) c write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) c write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) c write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) c call enerprint(energia(0)) write(iout,'(4hENER,i10,3f15.3,2i10,i5)') & i,energia(0),eini(i,iprot),entfac(i,iprot) write (iout,'(i5,18(1pe12.4))') ii, & (enetb(ii,j,iprot),j=1,n_ene) c write (iout,'(i5,16(1pe12.4))') i, c & (energia(j),j=1,n_ene),energia(0) c print *,"Processor",me,me1," computing energies",i #endif if (energia(0).ge.1.0d99) then write (iout,*) "FUNC1:CHUJ NASTAPIL in energy evaluation for", & " point",i,". Probably NaNs in some of the energy components." write (iout,*) "The components of the energy are:" call enerprint(energia(0)) write (iout,*) "Conformation",i,ii write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres) write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct) write (iout,'(8f10.4)') (vbld(k),k=2,nres) write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct) write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres) write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres) write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1) write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1) write (iout,*) "Calculation terminated at this point.", & " Check the database of conformations" call flush(iout) #ifdef MPI call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR) #endif stop "SEVERE error in energy calculation" endif else if (mod_side) then enetb(ii,1,iprot)=0.0d0 do j=1,ntyp do k=1,j enetb(ii,1,iprot)=enetb(ii,1,iprot)+ & ftune_eps(eps(j,k))*eneps(1,icant(j,k),ii,iprot)+ & eps(j,k)*eneps(2,icant(j,k),ii,iprot) enddo enddo endif etoti=0.0d0 do j=1,n_ene etoti=etoti+enetb(ii,j,iprot)*ww(j) enddo e_total(i,iprot)=etoti c write (iout,*) "FUNC1 natlike",iprot,natlike(iprot) if (natlike(iprot).gt.0) then call restore_ccoords(iprot,ii) call int_from_cart1(.false.) c call natsimil(i,ii,iprot,.false.) c write (iout,*) c & "natsimil",i,ii,(nu(k,ii,iprot),k=1,natlike(iprot)) endif #ifdef DEBUG write (iout,'(i5,18(1pe12.4))') i, & (enetb(ii,j,iprot),j=1,n_ene) write (iout,'(18(1pe12.4))') & (eneps(1,j,ii,iprot),j=201,210) write(iout,'(4hENER,i10,3f15.3)') & i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot) #endif endif return end