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" include "COMMON.DERIV" 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 #ifdef MPI double precision e_total_(maxstr_proc),eini_(maxstr_proc), & entfac_(maxstr_proc),rmstb_(maxstr_proc) #endif 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 lprn=.true. #ifdef DEBUG write (iout,*) "FUNC1: Init_Ene ",init_ene call flush(iout) write (iout,*) "ELOWEST at the beginning of FUC1" do iprot=1,nprot if (.not.maxlik(iprot)) cycle 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 c write (iout,*) "What",what," cutoffeval ",cutoffeval c call flush(iout) 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) c write (iout,*) "What",what," cutoffeval ",cutoffeval c call flush(iout) 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. #ifdef DEBUG write (iout,*) "Caling x2w" call flush(iout) #endif 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 if (.not.maxlik(iprot)) cycle #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," mod_ang",mod_ang #endif calc_grad=.false. 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) c Change AL 12/30/2017 if (.not. mod_other_params) & 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 .and. .not. mod_ang) & 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 c AL 12/30/17 - writing energies is not needed when whole energy is evaluated c 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 .and. .not. mod_ang) & 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 .or. mod_ang) & close(ientout) ENDIF #ifdef MPI c print *,"Processor",me,me1," finished computing energies" c call flush(iout) do i=1,scount(me1,iprot) e_total_(i)=e_total(indstart(me1,iprot)+i-1,iprot) eini_(i)=eini(indstart(me1,iprot)+i-1,iprot) entfac_(i)=entfac(indstart(me1,iprot)+i-1,iprot) rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot) enddo if (What.eq.2 .or. What.eq.3) then c call MPI_Gatherv(e_total(indstart(me1,iprot),iprot), c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot), c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, c & Comm1, IERROR) call MPI_Gatherv(e_total_(1), & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) c call MPI_Gatherv(entfac(indstart(me1,iprot),iprot), c & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot), c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, c & Comm1, IERROR) call MPI_Gatherv(entfac_(1), & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) c call MPI_Gatherv(eini(indstart(me1,iprot),iprot), c & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot), c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, cc & Comm1, IERROR) call MPI_Gatherv(eini_(1), & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot), & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, & Comm1, IERROR) c call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot), c & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot), c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master, c & Comm1, IERROR) call MPI_Gatherv(rmstb_(1), & 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 c call MPI_Allgatherv(e_total(indstart(me1,iprot),iprot), c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot), c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION, c & Comm1, IERROR) call MPI_Allgatherv(e_total_(1), & 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 C Change AL 12/30/2017 if (.not.mod_other_params) &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 if (.not.mod_other_params) &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 c write (iout,*) "Exit func1" 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 if (.not.maxlik(iprot)) cycle #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 if (.not.maxlik(iprot)) cycle 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 if (.not.maxlik(iprot)) cycle #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),g(max_paropt) 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" include "COMMON.TORSION" double precision ff common /sru/ ff integer i,ii,it,inat,j,k,l,kk,iprot,ib logical bs external ufparm double precision urparm,rdif integer uiparm double precision f,fpmf,chisquare_force,kara,viol,gnmr,gnmr1,aux double precision fforce logical lprn lprn = .true. #ifdef DEBUG write (iout,*) me,me1,"Calling TARGETFUN" write (iout,*) "n",n," x",(x(i),i=1,n) #endif call write_restart(n,x) nfl=nf call x2w(n,x) f=0.0d0 call func1(n,1,x) #ifdef NEWCORR if (mod_fourier(nloctyp).gt.0) then fpmf = rdif(n,x,g,1) else fpmf=0.0d0 endif #endif c write (iout,*) "Calling chisquare_force" c call flush(iout) fforce = chisquare_force(n,x,g,1) c write (iout,*) "Exit chisquare_force",fforce c call flush(iout) c Maximum likelihood contribution do iprot=1,nprot if (maxlik(iprot)) then do iT=1,nbeta(1,iprot) f = f + weilik(iT,iprot)*sumlik(iT,iprot) enddo endif enddo #ifdef DEBUG write (iout,*)"target_fun: sumlik" do iprot=1,nprot write (iout,*) (sumlik(ii,iprot),ii=1,nbeta(1,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,iprot), & " 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 f = f + wpmf*fpmf + fforce ff = f #ifdef DEBUG write (iout,*) "target_fun: f=",f,wpmf*fpmf,1.0d16*kara call flush(iout) #endif c#undef DEBUG 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),gforc(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" include "COMMON.TORSION" double precision ff common /sru/ ff #ifdef MPI include "COMMON.MPI" #endif double precision xi,delta /1.0d-6/ integer uiparm double precision urparm,rdif,chisquare_force 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),gang(maxang_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,nnene integer iti,itj,jj,icant integer l,ll logical lprn double precision gnmrprim,gnmr1prim double precision tcpu_ini,tcpu,chi0,chiplus,chiminus logical in_sc_pair_list external in_sc_pair_list integer ind_shield /25/ 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) #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then aux = rdif(n,x,gg,2) c write (iout,*) "After rdif" c call flush(iout) endif #endif chi0 = chisquare_force(n,x,gforc,2) 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 if (mod_ang) then do k=1,maxang_var gang(k)=0.0d0 enddo endif c Maximum likelihood gradient kk=0 do k=1,n_ene if (imask(k).gt.0 .and. k.ne.ind_shield) then kk=kk+1 endif enddo nnene=kk c write (iout,*) "kk",nn," nnene",nnene," n_ene",n_ene do iprot=1,nprot do ib=1,nbeta(1,iprot) kk=0 do k=1,n_ene if (imask(k).gt.0 .and. k.ne.ind_shield) 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 do k=1,nang_var gang(k)=gang(k)+weilik(ib,iprot)*sumlikang(k,ib,iprot) enddo enddo enddo #ifdef DEBUG write (iout,*) "gtor" do k=1,ntor_var write(iout,*) k,gtor(k) enddo write (iout,*) "nang_var",nang_var," gang" do k=1,nang_var write(iout,*) k,gang(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 .and. k.ne.ind_shield) 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 do k=1,nang_var gang(k)=gang(k)+zewnvtot(ib,iprot)*zvangtot(k,ib,iprot) enddo enddo enddo #ifdef DEBUG write (iout,*) "gtor" do k=1,ntor_var write(iout,*) k,gtor(k) enddo write (iout,*) "gang" do k=1,nang_var write(iout,*) k,gang(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 .and. k.ne.ind_shield) 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 do k=1,nang_var gang(k)=gang(k)+aux*rang(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 write (iout,*) "gang" do k=1,nang_var write(iout,*) k,gang(k) enddo #endif c Compute numerically the remaining part of the gradient c n_weight=0 c do i=1,n_ene c if (imask(i).gt.0 .and. i.ne.ind_shield) n_weight=n_weight+1 c enddo c ii=n_weight+17 ii=n_weight_PMF 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 if (mod_ang) then do i=1,nang_var ii=ii+1 g(ii)=gang(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) #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then aux = rdif(n,x,gg,0) c write (iout,*) "After rdif" c call flush(iout) endif #endif chi0 = chisquare_force(n,x,gforc,1) 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 c write (iout,*) "n_weight",n_weight," n",n c AL 5/12/19 numerical derivatives of force chisqaure in angles do i=nnene+1,n_weight 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) #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then aux = rdif(n,x,gg,0) #endif endif chiplus = chisquare_force(n,x,gforc,1) gforc(i)=-(chiplus-chi0)/delta x(i)=xi enddo c end nunumerical derivatives of forces in angles 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) #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then aux = rdif(n,x,gg,0) c write (iout,*) "After rdif" c call flush(iout) endif #endif chiplus = chisquare_force(n,x,gforc,1) c x(i)=xi-delta c call x2w(n,x) c call func1(n,1,x) c chiminus = chisquare_force(n,x,gforc,1) c write (iout,*) "chi0",chi0," chiplus",chiplus gforc(i)=-(chiplus-chi0)/delta c gforc(i)=-0.5d0*(chiplus-chiminus)/delta c Maximum likelihood do iprot=1,nprot do ib=1,nbeta(1,iprot) g(i)=g(i)+weilik(ib,iprot)*(sumlik(ib,iprot) & -sumlik0(ib,iprot))/delta c write (iout,*) "iprot",iprot," g",g(i) 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 NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then aux = rdif(n,x,gg,2) c write (iout,*) "After rdif" c call flush(iout) endif aux = chisquare_force(n,x,gforc,0) #endif #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 #ifdef NEWCORR c AL 2/10/2018 Add PMF gradient if (mod_fourier(nloctyp).gt.0) then do i=1,n g(i)=g(i)-wpmf*gg(i) enddo endif do i=1,n g(i)=g(i)-gforc(i) enddo #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), & enangavefj(0:3,maxang_var,2),enangavefjprim(0:3,maxang_var,2), & enangmixfj(0:3,maxang_var,2),enangmixsqfj(0:3,maxang_var,2) integer i,ii,j,k,kk,iprot,ib integer iti,itj,jj,icant integer l,ll logical lprn integer ind_shield /25/ 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 .and. k.ne.ind_shield) 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 if (mod_ang) then do k=1,nang_var ebisdecum=0.0d0 edecum=0.0d0 edecum1=enang_mix_ftot(k,ib,iprot) & -enangave_ftot(k,ib,iprot)*emean_ftot(ib,iprot) e2decum=enang_mixsq_ftot(k,ib,iprot) & -enangave_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, & " enangave_ftot",enangave_ftot(k,ib,iprot), & " enang_mix_ftot",enang_mix_ftot(k,ib,iprot), & " enang_mixsq_ftot",enang_mixsq_ftot(k,ib,iprot), & " emean_ftot",emean_ftot(ib,iprot) write (iout,*) ib,k," ebisdecum",ebisdecum, & " edecum",edecum," edecum1",edecum1," e2decum",e2decum, & " zvtortot",zvangtot(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 integer ind_shield /25/ 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 .and. k.ne.ind_shield) 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 if (mod_ang) then do l=1,natdim(inat,iprot) do k=1,nang_var rang(k,l)=-fac*(nuangave_f(k,l,ib,inat,iprot)- & nuave(l,inat,ib,iprot)*enangave_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" include "COMMON.LOCAL" include "COMMON.FFIELD" 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 integer ind_shield /25/ 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 & .or. mod_ang .or. mod_side_other .or. imask(ind_shield).gt.0) & ) 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 c#define DEBUG #ifdef DEBUG write (iout,*) "ww",ww(:n_ene) 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,e_total(i,iprot),eini(i,iprot),entfac(i,iprot) call enerprint(energia(0)) #endif c#undef DEBUG if (mod_side) then c 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 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 c AL 8/7/17: optimization of the bending potentials k=0 if (mod_ang) then do ij=0,nthetyp-1 if (mask_ang(ij).eq.0) cycle do jj=1,nbend_kcc_Tb(ij) k=k+1 enbend(k,ii,iprot)=ebend_temp_kcc(jj,ij)*wang c write (iout,*) "ij",ij," jj",jj," k",k," ebend_temp_kcc", c & ebend_temp_kcc(jj,ij) enddo enddo endif c#define DEBUG #ifdef DEBUG 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) call enerprint(energia(0)) c write(iout,'(4hENER,i10,3f15.3,2i10,i5)') c & i,energia(0),eini(i,iprot),entfac(i,iprot) c write (iout,'(i5,18(1pe12.4))') ii, c & (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 c#undef DEBUG 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 c#define DEBUG #ifdef DEBUG write (iout,*) "i",i," ww",ww write (iout,*) "enetb",(enetb(ii,j,iprot),j=1,n_ene) #endif c#undef DEBUG 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