1 subroutine func1(n,What,x)
4 include "DIMENSIONS.ZSCOPT"
12 include "COMMON.WEIGHTS"
13 include "COMMON.WEIGHTDER"
14 include "COMMON.CLASSES"
15 include "COMMON.ENERGIES"
16 include "COMMON.IOUNITS"
17 include "COMMON.VMCPAR"
18 include "COMMON.OPTIM"
19 include "COMMON.NAMES"
20 include "COMMON.INTERACT"
21 include "COMMON.TIME1"
22 include "COMMON.CHAIN"
23 include "COMMON.PROTFILES"
24 include "COMMON.COMPAR"
25 include "COMMON.PROTNAME"
26 include "COMMON.THERMAL"
27 include "COMMON.TORSION"
28 C Define local variables
29 integer i,ii,iii,ind,jj,j,k,l,ik,jk,iprot,ib,l1,l2,ll1,ll2,
32 integer ipass_conf,istart_conf,iend_conf
33 double precision energia(0:max_ene)
34 double precision etoti,sumpart,esquarei,emeani,elowesti,enepsjk,
35 & eave_pft(max_ene),emix_pft(max_ene),
36 & elowest_t(2,maxT),esquare_ft,
37 & elowest_aux(2,maxT),
38 & efree_t,emixsq_pft(max_ene),
39 & eneps_mixt(nntyp),eneps_meant(nntyp),
40 & enepsave_ft(nntyp),eneps_mix_ft(nntyp),
41 & eneps_mixsq_ft(nntyp),emean_ft,fac
43 double precision tcpu_ini,tcpu
51 c write (iout,*) "FUNC1: Init_Ene ",init_ene
55 write (iout,*) "ELOWEST at the beginning of FUC1"
57 do ibatch=1,natlike(iprot)+2
58 do ib=1,nbeta(ibatch,iprot)
59 write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
60 & " elowest",elowest(ib,ibatch,iprot)
66 call MPI_Bcast( What, 1, MPI_INTEGER, Master, Comm1, IERROR)
67 call MPI_Bcast( cutoffeval, 1, MPI_LOGICAL, Master, Comm1, IERROR)
68 call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION, Master, Comm1, IERROR)
71 write (iout,*) "Processor",me,me1," What",What
72 write (iout,*) "Processor",me,me1," x",(x(i),i=1,n)
77 write (iout,*) "x",(x(i),i=1,n)
80 C Convert variables into weights and other UNRES energy parameters.
82 C Calculate the total energies.
84 c write (iout,*) "Processor",me," nprot",nprot
85 c Revise the list of conformations
87 write (iout,*) "MAKE_LIST called"
89 call MPI_Bcast( enecut(1), nprot, MPI_DOUBLE_PRECISION,
90 & Master, Comm1, IERROR)
91 call make_list(.true.,.false.,n,x)
97 open(icsa_rbank,file="corr."//protname(iprot),
98 & access="direct",form="formatted",recl=25+10*natlike(iprot))
99 vormat='(2i6,2i3,e15.5,'
100 if (natlike(iprot).lt.10) then
101 write(vormat(16:),'(i1,a)') maxdimnat*natlike(iprot),'f7.4/)'
103 write(vormat(16:),'(i2,a)') maxdimnat*natlike(iprot),'f7.4/)'
105 write (iout,*) "format ",vormat
110 write (iout,*) "Processor",me," iprot",iprot,
111 & " indstart",indstart(me1,iprot)," indend",indend(me1,iprot),
112 & " init_ene",init_ene," mod_fourier",mod_fourier(nloctyp),
113 & " mod_elec",mod_elec," mod_scp",mod_scp," mod_tor",mod_tor,
114 & " mod_sccor",mod_sccor
116 call restore_molinfo(iprot)
119 IF (NCHUNK_CONF(IPROT).EQ.1) THEN
122 do i=indstart(me1,iprot),indend(me1,iprot)
124 do i=1,ntot_work(iprot)
126 c write (iout,*) "i",i," ii",ii," calling ene_eval"
127 call ene_eval(iprot,i,ii)
128 c write (iout,*) "After ene_eval: i",i," ii",ii
131 write (icsa_rbank,vormat,rec=i)
132 & i,list_conf(i,iprot),
133 & e_total(i,iprot),((nu(l,k,ii,iprot),l=1,natdim(k,iprot)),
134 & k=1,natlike(iprot))
141 if (.not.init_ene .and. mod_fourier(nloctyp).eq.0
142 & .and. .not. mod_elec .and. .not. mod_scp
143 & .and. .not. mod_tor .and. .not. mod_sccor) then
144 open (ientin,file=benefiles(iprot),status="old",
145 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
146 if (natlike(iprot).gt.0)
147 & open (icbase,file=bprotfiles(iprot),status="old",
148 & form="unformatted",access="direct",recl=lenrec(iprot))
150 open (icbase,file=bprotfiles(iprot),status="old",
151 & form="unformatted",access="direct",recl=lenrec(iprot))
152 c write (iout,*) "lenrec_ene",lenrec_ene(iprot)
154 open (ientout,file=benefiles(iprot),status="unknown",
155 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
159 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
161 iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
163 do i=1,ntot_work(iprot),maxstr_proc
165 iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
169 c Read the chunk of conformations off a DA scratchfile.
171 if (.not. init_ene .and. mod_fourier(nloctyp).eq.0
172 & .and. .not. mod_elec .and. .not. mod_scp
173 & .and. .not. mod_tor .and. .not. mod_sccor) then
175 write (iout,*) "func1: calling daread_ene",istart_conf,
178 call daread_ene(iprot,istart_conf,iend_conf)
179 if (natlike(iprot).gt.0) then
181 write (iout,*) "Calling daread_coords",istart_conf,
184 call daread_ccoords(iprot,istart_conf,iend_conf)
186 do iii=istart_conf,iend_conf
187 call ene_eval(iprot,iii,ii)
190 write (iout,*) "func1: calling dawrite_ene",istart_conf,
193 call dawrite_ene(iprot,istart_conf,iend_conf,ientin)
195 call daread_ccoords(iprot,istart_conf,iend_conf)
196 do iii=istart_conf,iend_conf
197 call ene_eval(iprot,iii,ii)
199 call dawrite_ene(iprot,istart_conf,iend_conf,ientout)
203 do iii=istart_conf,iend_conf
204 write (icsa_rbank,vormat,rec=iii)
205 & iii,list_conf(iii,iprot),
206 & e_total(iii,iprot),
207 & ((nu(l,k,iii-istart_conf+1,iprot),l=1,natdim(k,iprot)),
208 & k=1,natlike(iprot))
214 if (.not.init_ene .and. mod_fourier(nloctyp).eq.0
215 & .and. .not. mod_elec .and. .not. mod_scp
216 & .and. .not. mod_tor .and. .not. mod_sccor) then
218 if (natlike(iprot).gt.0) close(icbase)
222 if (init_ene .or. mod_fourier(nloctyp).gt.0 .or. mod_elec
223 & .or. mod_scp .or. mod_tor .or. mod_sccor) close(ientout)
227 c print *,"Processor",me,me1," finished computing energies"
230 if (What.eq.2 .or. What.eq.3) then
231 call MPI_Gatherv(e_total(indstart(me1,iprot),iprot),
232 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
233 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
235 call MPI_Gatherv(entfac(indstart(me1,iprot),iprot),
236 & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot),
237 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
239 call MPI_Gatherv(eini(indstart(me1,iprot),iprot),
240 & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot),
241 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
243 call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
244 & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
245 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
249 call MPI_Allgatherv(e_total(indstart(me1,iprot),iprot),
250 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
251 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
256 c Mean and lowest energies of structural classes
258 IF (NCHUNK_CONF(IPROT).EQ.1) THEN
261 do i=1,ntot_work(iprot)
265 do i=indstart(me1,iprot),indend(me1,iprot)
269 istart_conf=indstart(me1,iprot)
270 iend_conf=indend(me1,iprot)
272 do i=1,ntot_work(iprot)
276 iend_conf=ntot_work(iprot)
279 do i=1,ntot_work(iprot)
280 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
284 open (ientin,file=benefiles(iprot),status="old",
285 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
286 c call daread_ene(iprot,istart_conf,iend_conf)
287 call emin_search(iprot)
291 open (ientin,file=benefiles(iprot),status="old",
292 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
295 do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
296 iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
298 do istart_conf=1,ntot_work(iprot),maxstr_proc
299 iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
302 c Read the chunk of energies and derivatives off a DA scratchfile.
304 ipass_conf=ipass_conf+1
305 do i=1,ntot_work(iprot)
309 do i=istart_conf,iend_conf
314 write (iout,*) "ipass_conf",ipass_conf,
315 & " istart_conf",istart_conf," iend_conf",iend_conf
316 do i=1,ntot_work(iprot)
317 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
321 call daread_ene(iprot,istart_conf,iend_conf)
322 call emin_search(iprot)
331 c Complete the calculation of the lowest energies over all classes and
332 c distribute the values to all procs
335 write (iout,*) "Processor",me,me1," calling MPI_Reduce: 2"
336 do k=1,natlike(iprot)+2
337 write (iout,*) "iprot",iprot," k",k
338 do ib=1,nbeta(k,iprot)
339 write (iout,*) "ib",ib
340 write (iout,'(7hELOWEST,10e15.3)')
341 & elowest(ib,k,iprot)
346 do ibatch=1,natlike(iprot)+2
347 do ib=1,nbeta(ibatch,iprot)
348 elowest_aux(1,ib)=elowest(ib,ibatch,iprot)
349 elowest_aux(2,ib)=ind_lowest(ib,ibatch,iprot)
351 call MPI_Allreduce(elowest_aux(1,1),elowest_t(1,1),
352 & nbeta(ibatch,iprot),
353 & MPI_2DOUBLE_PRECISION, MPI_MINLOC, Comm1, IERROR)
355 do ib=1,nbeta(ibatch,iprot)
356 write (iout,*) "beta=",betaT(ib,ibatch,iprot)
357 write (iout,'(9helowest_t,10f15.3)')
358 & elowest_t(1,ib),elowest_t(2,ib)
360 write (iout,*) "Processor",me,me1," finished MPI_Reduce: 2"
362 do ib=1,nbeta(ibatch,iprot)
363 c write (iout,*) "IB",ib," ibatch",ibatch," iprot",iprot
365 elowest(ib,ibatch,iprot)=elowest_t(1,ib)
366 c write (iout,*) "elowest",ib,ibatch,iprot,
367 c & elowest(ib,ibatch,iprot)
369 ind_lowest(ib,ibatch,iprot)=elowest_t(2,ib)
377 c write (iout,*) "After averages, calling Thermal"
379 if (What.eq.2) call thermal(iprot)
380 c write (iout,*) "After thermal"
384 if (me1.eq.Master) then
388 write (iout,*) "iprot",iprot," elowest",
389 & ((elowest(ib,k,iprot),ib=1,nbeta(k,iprot)),k=1,natlike(iprot)+2)
392 do ib=1,nbeta(2,iprot)
393 fac=betaT(ib,2,iprot)
394 c 6/15/05 AL Heat capacity
395 zvtot(ib,iprot)=(esquare_ftot(ib,iprot)-
396 & emean_ftot(ib,iprot)**2)*fac*fac*Rgas
397 & -ebis_ftot(ib,iprot)+heatbase(iprot)
399 write (iout,*) "iprot",iprot," ib",ib,
400 & " emeantot",emean_ftot(ib,iprot),
401 & " esquaretot",esquare_ftot(ib,iprot),
402 & " ebistot",ebis_ftot(ib,iprot),
403 & " zvtot",zvtot(ib,iprot)
406 write (iout,*) "beta",
408 & " efree",efree(ib,2,iprot),
409 & " emean",emean_ftot(ib,iprot),
410 & " variance",esquare_ftot(ib,iprot)
413 c 4/13/04 AL Correlation coefficients between energy and nativelikeness
414 c in the respective classes.
416 write (iout,*) "iprot",iprot," ib",ib,
418 write (iout,*) "emean",emean_ftot(ib,iprot),
419 & " esquare",esquare_ftot(ib,iprot)
425 & "Averages of nativelikeness and energy-nativelikeness",
426 & " correlation coefficients:"
428 do i=1,natlike(iprot)
429 do ib=1,nbeta(i+1,iprot)
430 write (iout,*) "protein",iprot," property",i,
431 & " temperature",ib," natdim",natdim(i,iprot)," value(s)",
432 & (nuave(j,ib,i,iprot),j=1,natdim(i,iprot))
443 write (iout,*) "Calling CUTOFF_VIOLATION"
446 c write (iout,*) "CUTOFFEVAL",cutoffeval
448 if (cutoffeval) call cutoff_violation
449 t_func = t_func + tcpu() - tcpu_ini
451 write (iout,*) "Exitting FUNC1"
456 c-------------------------------------------------------------------------------
457 subroutine targetfun(n,x,nf,f,uiparm,ufparm)
460 include "DIMENSIONS.ZSCOPT"
463 integer IERROR,Errcode,status(MPI_STATUS_SIZE)
467 double precision x(n)
468 include "COMMON.WEIGHTS"
469 include "COMMON.IOUNITS"
470 include "COMMON.DERIV"
471 include "COMMON.VMCPAR"
472 include "COMMON.CLASSES"
474 include "COMMON.ENERGIES"
475 include "COMMON.WEIGHTDER"
476 include "COMMON.OPTIM"
477 include "COMMON.TIME1"
478 include "COMMON.COMPAR"
481 integer i,ii,it,inat,j,k,l,kk,iprot,ib
484 double precision urparm
486 double precision f,kara,viol,gnmr,gnmr1,aux
490 write (iout,*) me,me1,"Calling TARGETFUN"
492 call write_restart(n,x)
498 c Maximum likelihood contribution
500 do iT=1,nbeta(1,iprot)
501 f = f + weilik(iT,iprot)*sumlik(iT,iprot)
505 write (iout,*)"target_fun: sumlik"
507 write (iout,*) (sumlik(ii,iprot),ii=1,nbeta_l(iprot))
509 write (iout,*) "f before Cv =",f
513 c 6/15/05 AL Contribution from Cv
514 do ib=1,nbeta(2,iprot)
515 f=f+weiCv(ib,iprot)*(zvtot(ib,iprot)-target_cv(ib,iprot))**2
517 c write (iout,*) "target_fun from CV: f=",f
519 c write (iout,*) "target_fun: f=",f
520 c 8/2/2013 Contribution from conformation-dependent averages
522 do inat=1,natlike(iprot)
523 do ib=1,nbeta(inat+2,iprot)
524 do k=1,natdim(inat,iprot)
525 f = f+weinat(k,ib,inat,iprot)*
526 & (nuave(k,ib,inat,iprot)-nuexp(k,ib,inat,iprot))**2
528 write (iout,*) "iprot",iprot," inat",inat,
529 & " ib",ib," k=",k," nuave",nuave(k,ib,inat,iproti),
530 & " nuexp",nuexp(k,ib,inat,iprot),
531 & " weight",weinat(k,ib,inat,iprot)
537 c write (iout,*) "target_fun after adding nativelikeness: f=",f
539 c add constraints on weights
542 kara=kara+gnmr1(x(i),x_low(i),x_up(i))
547 c write (iout,*) "target_fun: f=",f,1.0d16*kara
551 c-----------------------------------------------------------------------------
552 subroutine targetgrad(n,x,nf,g,uiparm,rdum,ufparm)
555 include "DIMENSIONS.ZSCOPT"
556 include "COMMON.TIME1"
559 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
562 double precision x(n),
563 & g(n),gg(max_paropt)
564 include "COMMON.WEIGHTS"
565 include "COMMON.ENERGIES"
566 include "COMMON.CLASSES"
567 include "COMMON.IOUNITS"
568 include "COMMON.DERIV"
569 include "COMMON.COMPAR"
570 include "COMMON.WEIGHTDER"
571 include "COMMON.OPTIM"
572 include "COMMON.VMCPAR"
574 include "COMMON.NAMES"
581 double precision xi,delta /1.0d-6/
583 double precision urparm
585 double precision zewn(MaxT,maxprot),aux,zscder
586 double precision zewnvtot(maxT,maxprot)
587 double precision sumlik0(maxT,maxprot),zv0tot(maxT,maxprot),
588 & nuave0(maxdimnat,maxT,maxnatlike,maxprot)
589 double precision geps(nntyp),gtor(maxtor_var)
590 double precision rdum
592 integer ll1,ll2,nl1,nl2,nn
593 integer i,ii,ind,inat,j,k,kk,iprot,ib,n_weight
594 integer iti,itj,jj,icant
597 double precision gnmrprim,gnmr1prim
598 double precision tcpu_ini,tcpu
599 logical in_sc_pair_list
600 external in_sc_pair_list
606 write (iout,*) "Processor",me,me1," calling TARGETGRAD"
611 write (iout,*) "Calling TARGETGRAD"
616 write (iout,*) "x",(x(i),i=1,n)
623 do ib=1,nbeta(2,iprot)
624 c 6/15/05 AL Contributions from Cv
625 zewnvtot(ib,iprot)=2*weiCv(ib,iprot)*
626 & (zvtot(ib,iprot)-target_cv(ib,iprot))
638 if (mod_tor .or. mod_sccor) then
643 c Maximum likelihood gradient
645 do ib=1,nbeta(1,iprot)
648 if (imask(k).gt.0) then
650 g(kk)=g(kk)+weilik(ib,iprot)*sumlikder(kk,ib,iprot)
654 geps(k)=geps(k)+weilik(ib,iprot)*sumlikeps(k,ib,iprot)
657 gtor(k)=gtor(k)+weilik(ib,iprot)*sumliktor(k,ib,iprot)
662 write (iout,*) "gtor"
664 write(iout,*) k,gtor(k)
668 c Heat-capacity gradient
670 do ib=1,nbeta(2,iprot)
673 if (imask(k).gt.0) then
675 g(kk)=g(kk)+zvdertot(kk,ib,iprot)*zewnvtot(ib,iprot)
679 geps(k)=geps(k)+zewnvtot(ib,iprot)*zvepstot(k,ib,iprot)
682 gtor(k)=gtor(k)+zewnvtot(ib,iprot)*zvtortot(k,ib,iprot)
687 write (iout,*) "gtor"
689 write(iout,*) k,gtor(k)
691 write (iout,*) "grad: g",(g(i),i=1,n)
694 c Derivatives of conformational-dependent properties
696 do inat=1,natlike(iprot)
697 do ib=1,nbeta(inat+2,iprot)
698 call propderiv(ib,inat,iprot)
699 do l=1,natdim(inat,iprot)
700 aux=2*weinat(l,ib,inat,iprot)*
701 & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
703 write (iout,*) "iprot",iprot," ib",ib,
704 & " inat",inat," l",l,
705 & " nuave",nuave(l,ib,inat,iprot),
706 & " weight",weinat(ib,inat,iprot),
712 if (imask(k).gt.0) then
714 g(kk)=g(kk)+aux*rder(kk,l)
718 geps(k)=geps(k)+aux*reps(k,l)
721 gtor(k)=gtor(k)+aux*rtor(k,l)
728 write (iout,*) "gtor"
730 write(iout,*) k,gtor(k)
733 c Compute numerically the remaining part of the gradient
736 if (imask(i).gt.0) n_weight=n_weight+1
739 c Sum up the gradient in SC epsilons, if needed
744 g(ii)=g(ii)+geps(icant(iti,iti))
746 if (.not. in_sc_pair_list(iti,j)) then
747 g(ii)=g(ii)+0.5d0*geps(icant(iti,j))
755 g(ii)=geps(icant(iti,itj))
758 c write (iout,*) "n_ene",n_ene
759 if (mod_tor .or. mod_sccor) then
767 write (iout,*) "n",n," n_weight",n_weight
768 write (iout,*) "After SC: grad:"
770 write (iout,*) i,g(i)
776 do ib=1,nbeta(1,iprot)
777 sumlik0(ib,iprot)=sumlik(ib,iprot)
781 do ib=1,nbeta(2,iprot)
782 zv0tot(ib,iprot)=zvtot(ib,iprot)
786 do i=1,natlike(iprot)
787 do ib=1,nbeta(i+2,iprot)
788 do k=1,natdim(i,iprot)
789 nuave0(k,ib,i,iprot)=nuave(k,ib,i,iprot)
796 if (nbeta(2,iprot).gt.0) nn=nn-1
800 if (nbeta(2,iprot).gt.0) then
802 do ib=1,nbeta(2,iprot)
803 g(ii)=g(ii)+zewnvtot(ib,iprot)
808 write (iout,*) "n",n," n_weight",n_weight
809 write (iout,*) "After heat grad:"
811 write (iout,*) i,g(i)
815 c write (iout,*) "n_weight",n_weight," n",n
819 c write (iout,*) "i",i," x",x(i),xi
824 do j=1,nbeta(1,iprot)
825 g(i)=g(i)+weilik(ib,iprot)*(sumlik(j,iprot)
826 & -sumlik0(j,iprot))/delta
830 write (iout,*) "init numerical derivatives"
833 c Contribution from Cv
835 do ib=1,nbeta(2,iprot)
836 g(i)=g(i)+zewnvtot(ib,iprot)*
837 & (zvtot(ib,iprot)-zv0tot(ib,iprot))/delta
840 c Contribution from average nativelikenss and correlation coefficients
841 c between energy and nativelikeness
843 do ii=1,natlike(iprot)
844 do ib=1,nbeta(ii+2,iprot)
845 do k=1,natdim(i,iprot)
846 aux=2*weinat(k,ib,inat,iprot)*
847 & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
848 g(i)=g(i)+aux*(nuave(k,ib,ii,iprot)-
849 & nuave0(k,ib,ii,iprot))/delta
857 write (iout,*) "n",n," n_weight",n_weight
858 write (iout,*) "After num grad:"
860 write (iout,*) k,g(k)
866 write (iout,*) "finish numerical derivatives"
870 c 1/10/2007 AL Correction: ENDIF...IF insterted to prevent unmatched
871 c SEND is I_LOCAL_CHECK.GT.0
874 write (iout,*) "mode",mode
877 c transform the gradient to the x space
879 g(i)=g(i)*(x_up(i)-x_low(i))/(pi*(1.0d0+(x(i)
880 & -0.5d0*(x_up(i)+x_low(i)))**2))
883 c add constraints on weights
884 c write (iout,*) "Gradient in constraints"
886 g(i)=g(i)+1.0d16*gnmr1prim(x(i),x_low(i),x_up(i))
888 write (iout,*) i,x(i),x_low(i),x_up(i),
889 & gnmr1prim(x(i),x_low(i),x_up(i))
894 write (iout,*) "grad"
896 write (iout,*) i,g(i)
899 t_grad = t_grad + tcpu() - tcpu_ini
901 write (iout,*) "Exit targetgrad"
906 c------------------------------------------------------------------------------
907 subroutine matout(ndim,n,a)
908 implicit real*8 (a-h,o-z)
909 include "DIMENSIONS.ZSCOPT"
910 include "COMMON.IOUNITS"
911 double precision a(ndim)
912 character*6 line6 / '------' /
913 character*12 line12 / '------------' /
917 write (iout,1000) (k,k=i,nlim)
918 write (iout,1020) line6,(line12,k=i,nlim)
919 1000 format (/7x,6(i7,5x))
920 1020 format (a6,6a12)
925 3 b(k-i+1)=a(icant(i,j))
926 write (iout,1010) j,(b(k),k=1,jlim1)
929 1010 format (i3,3x,6(1pe11.4,1x))
932 c------------------------------------------------------------------------------
936 include "DIMENSIONS.ZSCOPT"
939 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
942 include "COMMON.WEIGHTS"
943 include "COMMON.ENERGIES"
944 include "COMMON.IOUNITS"
945 include "COMMON.DERIV"
946 include "COMMON.WEIGHTDER"
947 include "COMMON.VMCPAR"
949 include "COMMON.NAMES"
951 include "COMMON.CLASSES"
952 include "COMMON.THERMAL"
956 include "COMMON.TIME1"
957 double precision aux,zscder
958 double precision fac,fac2
959 double precision efj(2),emj(2),varj(2),edecum,edecum1,ebisdecum,
961 double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
962 & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
963 & enepsmixsqfj(nntyp,2),eavpfjprim(max_ene,2),
964 & entoravefj(0:3,maxtor_var,2),entoravefjprim(0:3,maxtor_var,2),
965 & entormixfj(0:3,maxtor_var,2),entormixsqfj(0:3,maxtor_var,2)
966 integer i,ii,j,k,kk,iprot,ib
967 integer iti,itj,jj,icant
973 do ib=1,nbeta(2,iprot)
974 fac = betaT(ib,2,iprot)
978 if (imask(k).gt.0) then
980 ebisdecum=emix_pfbistot(k,ib,iprot)-
981 & eave_pftot(k,ib,iprot)
982 & *ebis_ftot(ib,iprot)
983 edecum=emix_pfprimtot(k,ib,iprot)
984 & -eave_pfprimtot(k,ib,iprot)
985 & *emean_ftot(ib,iprot)
986 edecum1=emix_pftot(k,ib,iprot)
987 & -eave_pftot(k,ib,iprot)
988 & *emean_ftot(ib,iprot)
989 e2decum=emixsq_pftot(k,ib,iprot)
990 & -eave_pftot(k,ib,iprot)
991 & *esquare_ftot(ib,iprot)
992 zvdertot(kk,ib,iprot)=Rgas*fac2*(2*edecum
993 & -fac*(e2decum-2*emean_ftot(ib,iprot)
995 & -eave_pfbistot(k,ib,iprot)
998 write (iout,*) "k=",k,
999 & "eave_pftot",eave_pftot(k,ib,iprot),
1000 & " emix_pftot",emix_pftot(k,ib,iprot),
1001 & " emix_pfprimtot",emix_pfprimtot(k,ib,iprot),
1002 & " emix_pfbistot",emix_pfbistot(k,ib,iprot),
1003 & " eave_pfprimtot",eave_pfprimtot(k,ib,iprot),
1004 & " eave_pfbistot",eave_pfbistot(k,ib,iprot),
1005 & " emixsq_pftot",emixsq_pftot(k,ib,iprot),
1006 & " ebis_ftot",ebis_ftot(ib,iprot),
1007 & " emean_ftot",emean_ftot(ib,iprot)
1008 write (iout,*) "k=",k," edecum",edecum,
1009 & " edecum1",edecum1,
1010 & " e2decum",e2decum," ebisdecum",ebisdecum,
1012 & zvdertot(kk,ib,iprot)
1018 ebisdecum=eneps_mix_fbistot(k,ib,iprot)-
1019 & enepsave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1020 edecum=eneps_mix_ftot(k,ib,iprot)
1021 & -enepsave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1022 e2decum=eneps_mixsq_ftot(k,ib,iprot)
1023 & -enepsave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1024 zvepstot(k,ib,iprot)=ww(1)*(Rgas*fac2*(2*edecum
1025 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1026 & *edecum))+fac*ebisdecum)
1028 write (iout,*) ib,k,"ebisdecum",ebisdecum,
1029 & " edecum",edecum," e2decum",e2decum," zvepstot",
1030 & zvepstot(k,ibb,iprot)
1032 & " enepsave_ftot",enepsave_ftot(k,ib,iprot),
1033 & " eneps_mix_ftot",eneps_mix_ftot(k,ib,iprot),
1034 & " eneps_mix_fbistot",eneps_mix_fbistot(k,ib,iprot),
1035 & " eneps_mixsq_ftot",eneps_mixsq_ftot(k,ib,iprot),
1036 & " emean_ftot",emean_ftot(ib,iprot),
1037 & " ebis_ftot",ebis_ftot(ib,iprot)
1041 if (mod_tor .or. mod_sccor) then
1043 ebisdecum=entor_mix_fbistot(k,ib,iprot)-
1044 & entorave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1045 edecum=entor_mix_fprimtot(k,ib,iprot)
1046 & -entorave_fprimtot(k,ib,iprot)*emean_ftot(ib,iprot)
1047 edecum1=entor_mix_ftot(k,ib,iprot)
1048 & -entorave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1049 e2decum=entor_mixsq_ftot(k,ib,iprot)
1050 & -entorave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1051 zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1052 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1053 & *edecum1))-entorave_fbistot(k,ib,iprot)
1057 & " entorave_ftot",entorave_ftot(k,ib,iprot),
1058 & " entor_mix_ftot",entor_mix_ftot(k,ib,iprot),
1059 & " entor_mix_fprimtot",entor_mix_fprimtot(k,ib,iprot),
1060 & " entor_mix_fbistot",entor_mix_fbistot(k,ib,iprot),
1061 & " entorave_fprimtot",entorave_fprimtot(k,ib,iprot),
1062 & " entorave_fbistot",entorave_fbistot(k,ib,iprot),
1063 & " entor_mixsq_ftot",entor_mixsq_ftot(k,ib,iprot),
1064 & " emean_ftot",emean_ftot(ib,iprot),
1065 & " ebis_ftot",ebis_ftot(ib,iprot)
1066 write (iout,*) ib,k," ebisdecum",ebisdecum,
1067 & " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1068 & " zvtortot",zvtortot(k,ibb,iprot)
1073 write (iout,*) "iprot",iprot,
1075 write (iout,*) "zvder",
1076 & (zvdertot(k,ibb,iprot),k=1,kk)
1082 c------------------------------------------------------------------------------
1083 subroutine propderiv(ib,inat,iprot)
1084 c Compute the derivatives of conformation-dependent averages in force-field
1087 include "DIMENSIONS"
1088 include "DIMENSIONS.ZSCOPT"
1091 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1094 include "COMMON.WEIGHTS"
1095 include "COMMON.ENERGIES"
1096 include "COMMON.CLASSES"
1097 include "COMMON.IOUNITS"
1098 include "COMMON.DERIV"
1099 include "COMMON.WEIGHTDER"
1100 include "COMMON.VMCPAR"
1101 include "COMMON.GEO"
1102 include "COMMON.OPTIM"
1103 include "COMMON.NAMES"
1104 include "COMMON.COMPAR"
1106 include "COMMON.MPI"
1108 include "COMMON.TIME1"
1109 include "COMMON.VAR"
1110 double precision aux,zscder
1111 double precision fac
1112 double precision efj(2),emj(2),varj(2)
1113 double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1114 & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1115 & enepsmixsqfj(nntyp,2)
1117 integer i,ii,ind,j,k,kk,iprot,ib,ibb,l1
1118 integer iti,itj,jj,icant
1121 c write (iout,*) "Entered PROPDERIV"
1123 fac = betaT(ib,inat+2,iprot)
1124 c write (iout,*) "derivatives: iprot",iprot,
1126 do l=1,natdim(inat,iprot)
1129 if (imask(k).gt.0) then
1131 rder(kk,l)=-fac*(nu_pf(k,l,ib,inat,iprot)
1132 & -nuave(l,ib,inat,iprot)*
1133 & eave_nat_pftot(k,ib,inat,iprot))
1138 do l=1,natdim(inat,iprot)
1140 reps(k,l)=-ww(1)*fac*(nuepsave_f(k,l,ib,inat,iprot)-
1141 & nuave(l,ib,inat,iprot)*enepsave_nat_ftot(k,ib,inat,iprot))
1145 if (mod_tor .or. mod_sccor) then
1146 do l=1,natdim(inat,iprot)
1148 rtor(k,l)=-fac*(nutorave_f(k,l,ib,inat,iprot)-
1149 & nuave(l,inat,ib,iprot)*entorave_nat_ftot(k,ib,inat,iprot))
1155 c-------------------------------------------------------------------------------
1156 subroutine ene_eval(iprot,i,ii)
1158 include "DIMENSIONS"
1159 include "DIMENSIONS.ZSCOPT"
1162 integer IERROR,ErrCode
1163 include "COMMON.MPI"
1165 include "COMMON.CONTROL"
1166 include "COMMON.COMPAR"
1167 include "COMMON.WEIGHTS"
1168 include "COMMON.WEIGHTDER"
1169 include "COMMON.CLASSES"
1170 include "COMMON.ENERGIES"
1171 include "COMMON.IOUNITS"
1172 include "COMMON.VMCPAR"
1173 include "COMMON.NAMES"
1174 include "COMMON.INTERACT"
1175 include "COMMON.TIME1"
1176 include "COMMON.CHAIN"
1177 include "COMMON.VAR"
1178 include "COMMON.GEO"
1179 include "COMMON.SCCOR"
1180 include "COMMON.TORSION"
1181 C Define local variables
1182 integer i,ii,ij,jj,j,k,ik,jk,l,ll,lll,iprot,iblock
1183 double precision energia(0:max_ene)
1184 double precision etoti,enepsjk
1185 double precision ftune_eps
1193 write (iout,*) "ene_eval: iprot",iprot," i",i," ii",ii
1195 if (init_ene .or. (mod_fourier(nloctyp).gt.0
1196 & .or. mod_elec .or. mod_scp .or. mod_tor .or. mod_sccor) ) then
1198 write (iout,*) "FUNC1: doing UNRES"
1199 write (iout,*) "Protein",iprot," i",i," ii",ii," nres",nres
1202 call restore_ccoords(iprot,ii)
1203 call int_from_cart1(.false.)
1205 write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1206 write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1207 write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1208 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1209 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1210 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1211 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1212 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1217 call etotal(energia(0))
1218 c if (natlike(iprot).gt.0)
1219 c & call natsimil(i,ii,iprot,.false.)
1220 e_total(i,iprot)=energia(0)
1222 enetb(ii,j,iprot)=energia(j)
1225 write (iout,*) "FUN nntyp",nntyp," eneps",eneps_temp(1,45)
1228 eneps(k,j,ii,iprot)=eneps_temp(k,j)
1232 write (iout,'(2i5,18(1pe12.4))') i,ii,
1233 & (enetb(ii,j,iprot),j=1,n_ene)
1234 write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1235 & i,energia(0),eini(i,iprot),entfac(i,iprot)
1238 c write (iout,*) "Matrix eneps"
1241 c write (iout,*) j,k,icant(j,k),
1242 c & eneps_temp(1,icant(j,k)),
1243 c & eneps_temp(2,icant(j,k))
1249 if (mod_tor .or. mod_sccor) then
1250 c write (iout,*) "etor",enetb(ii,13,iprot),enetb(ii,19,iprot)
1253 do ij=-ntortyp+1,ntortyp-1
1254 do jj=-ntortyp+1,ntortyp-1
1255 if (mask_tor(0,jj,ij,iblock).eq.1) then
1257 entor(k,ii,iprot)=etor_temp(0,0,jj,ij,iblock)
1258 c write (iout,*) "etor_temp",restyp(itortyp(j)),
1259 c & restyp(itortyp(i)),
1260 c & k,etor_temp(ji,ii)
1261 else if (mask_tor(0,jj,ij,iblock).gt.1) then
1262 if (tor_mode.eq.0) then
1263 do l=1,nterm(jj,ij,iblock)
1265 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1267 if (mask_tor(0,jj,ij,iblock).gt.2) then
1268 do l=nterm(jj,ij,iblock)+1,2*nterm(jj,ij,iblock)
1270 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1273 else if (tor_mode.eq.1) then
1275 do l=1,nterm_kcc(jj,ij)
1276 do ll=1,nterm_kcc_Tb(jj,ij)
1277 do lll=1,nterm_kcc_Tb(jj,ij)
1279 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1283 if (mask_tor(0,jj,ij,1).gt.2) then
1284 do l=nterm_kcc(jj,ij)+1,2*nterm_kcc(jj,ij)
1285 do ll=1,nterm_kcc_Tb(jj,ij)
1286 do lll=1,nterm_kcc_Tb(jj,ij)
1288 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1299 do ij=-nsccortyp+1,nsccortyp-1
1300 do jj=-nsccortyp+1,nsccortyp-1
1302 if (mask_tor(ll,jj,ij,1).eq.1) then
1304 entor(k,ii,iprot)=etor_temp(0,ll,jj,ij,1)
1305 c write (iout,*) "etor_temp",restyp(isccortyp(jj)),
1306 c & restyp(isccortyp(ij)),
1307 c & k,etor_temp(jj,ij)
1308 else if (mask_tor(ll,jj,ij,1).gt.1) then
1309 do l=1,nterm_sccor(jj,ij)
1311 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1313 do l=nterm_sccor(jj,ij)+1,2*nterm_sccor(jj,ij)
1315 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1323 c write (iout,*) "Conformation",i
1324 c write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1325 c write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1326 c write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1327 c write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1328 c write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1329 c write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1330 c write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1331 c write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1332 c call enerprint(energia(0))
1333 write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1334 & i,energia(0),eini(i,iprot),entfac(i,iprot)
1335 write (iout,'(i5,18(1pe12.4))') ii,
1336 & (enetb(ii,j,iprot),j=1,n_ene)
1337 c write (iout,'(i5,16(1pe12.4))') i,
1338 c & (energia(j),j=1,n_ene),energia(0)
1339 c print *,"Processor",me,me1," computing energies",i
1341 if (energia(0).ge.1.0d99) then
1342 write (iout,*) "FUNC1:CHUJ NASTAPIL in energy evaluation for",
1343 & " point",i,". Probably NaNs in some of the energy components."
1344 write (iout,*) "The components of the energy are:"
1345 call enerprint(energia(0))
1346 write (iout,*) "Conformation",i,ii
1347 write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1348 write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1349 write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1350 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1351 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1352 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1353 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1354 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1355 write (iout,*) "Calculation terminated at this point.",
1356 & " Check the database of conformations"
1359 call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
1361 stop "SEVERE error in energy calculation"
1365 enetb(ii,1,iprot)=0.0d0
1368 enetb(ii,1,iprot)=enetb(ii,1,iprot)+
1369 & ftune_eps(eps(j,k))*eneps(1,icant(j,k),ii,iprot)+
1370 & eps(j,k)*eneps(2,icant(j,k),ii,iprot)
1376 etoti=etoti+enetb(ii,j,iprot)*ww(j)
1378 e_total(i,iprot)=etoti
1379 c write (iout,*) "FUNC1 natlike",iprot,natlike(iprot)
1380 if (natlike(iprot).gt.0) then
1381 call restore_ccoords(iprot,ii)
1382 call int_from_cart1(.false.)
1383 c call natsimil(i,ii,iprot,.false.)
1385 c & "natsimil",i,ii,(nu(k,ii,iprot),k=1,natlike(iprot))
1388 write (iout,'(i5,18(1pe12.4))') i,
1389 & (enetb(ii,j,iprot),j=1,n_ene)
1390 write (iout,'(18(1pe12.4))')
1391 & (eneps(1,j,ii,iprot),j=201,210)
1392 write(iout,'(4hENER,i10,3f15.3)')
1393 & i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)