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 e_total_(maxstr_proc),eini_(maxstr_proc),
44 & entfac_(maxstr_proc),rmstb_(maxstr_proc)
47 double precision tcpu_ini,tcpu
55 c write (iout,*) "FUNC1: Init_Ene ",init_ene
59 write (iout,*) "ELOWEST at the beginning of FUC1"
61 do ibatch=1,natlike(iprot)+2
62 do ib=1,nbeta(ibatch,iprot)
63 write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
64 & " elowest",elowest(ib,ibatch,iprot)
70 call MPI_Bcast( What, 1, MPI_INTEGER, Master, Comm1, IERROR)
71 call MPI_Bcast( cutoffeval, 1, MPI_LOGICAL, Master, Comm1, IERROR)
72 call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION, Master, Comm1, IERROR)
75 write (iout,*) "Processor",me,me1," What",What
76 write (iout,*) "Processor",me,me1," x",(x(i),i=1,n)
81 write (iout,*) "x",(x(i),i=1,n)
84 C Convert variables into weights and other UNRES energy parameters.
86 C Calculate the total energies.
88 c write (iout,*) "Processor",me," nprot",nprot
89 c Revise the list of conformations
91 if (me.eq.Master) then
92 write (iout,*) "MAKE_LIST called"
95 call MPI_Bcast( enecut(1), nprot, MPI_DOUBLE_PRECISION,
96 & Master, Comm1, IERROR)
97 call make_list(.true.,.false.,n,x)
103 open(icsa_rbank,file="corr."//protname(iprot),
104 & access="direct",form="formatted",recl=25+10*natlike(iprot))
105 vormat='(2i6,2i3,e15.5,'
106 if (natlike(iprot).lt.10) then
107 write(vormat(16:),'(i1,a)') maxdimnat*natlike(iprot),'f7.4/)'
109 write(vormat(16:),'(i2,a)') maxdimnat*natlike(iprot),'f7.4/)'
111 c write (iout,*) "format ",vormat
116 write (iout,*) "Processor",me," iprot",iprot,
117 & " indstart",indstart(me1,iprot)," indend",indend(me1,iprot),
118 & " init_ene",init_ene," mod_fourier",mod_fourier(nloctyp),
119 & " mod_elec",mod_elec," mod_scp",mod_scp," mod_tor",mod_tor,
120 & " mod_sccor",mod_sccor
122 call restore_molinfo(iprot)
125 IF (NCHUNK_CONF(IPROT).EQ.1) THEN
128 do i=indstart(me1,iprot),indend(me1,iprot)
130 do i=1,ntot_work(iprot)
132 c write (iout,*) "i",i," ii",ii," calling ene_eval"
133 call ene_eval(iprot,i,ii)
134 c write (iout,*) "After ene_eval: i",i," ii",ii
137 write (icsa_rbank,vormat,rec=i)
138 & i,list_conf(i,iprot),
139 & e_total(i,iprot),((nu(l,k,ii,iprot),l=1,natdim(k,iprot)),
140 & k=1,natlike(iprot))
147 if (.not.init_ene .and. mod_fourier(nloctyp).eq.0
148 & .and. .not. mod_elec .and. .not. mod_scp
149 & .and. .not. mod_tor .and. .not. mod_sccor) then
150 open (ientin,file=benefiles(iprot),status="old",
151 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
152 if (natlike(iprot).gt.0)
153 & open (icbase,file=bprotfiles(iprot),status="old",
154 & form="unformatted",access="direct",recl=lenrec(iprot))
156 open (icbase,file=bprotfiles(iprot),status="old",
157 & form="unformatted",access="direct",recl=lenrec(iprot))
158 c write (iout,*) "lenrec_ene",lenrec_ene(iprot)
160 open (ientout,file=benefiles(iprot),status="unknown",
161 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
165 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
167 iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
169 do i=1,ntot_work(iprot),maxstr_proc
171 iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
175 c Read the chunk of conformations off a DA scratchfile.
177 if (.not. init_ene .and. mod_fourier(nloctyp).eq.0
178 & .and. .not. mod_elec .and. .not. mod_scp
179 & .and. .not. mod_tor .and. .not. mod_sccor) then
181 write (iout,*) "func1: calling daread_ene",istart_conf,
184 call daread_ene(iprot,istart_conf,iend_conf)
185 if (natlike(iprot).gt.0) then
187 write (iout,*) "Calling daread_coords",istart_conf,
190 call daread_ccoords(iprot,istart_conf,iend_conf)
192 do iii=istart_conf,iend_conf
193 call ene_eval(iprot,iii,ii)
196 write (iout,*) "func1: calling dawrite_ene",istart_conf,
199 call dawrite_ene(iprot,istart_conf,iend_conf,ientin)
201 call daread_ccoords(iprot,istart_conf,iend_conf)
202 do iii=istart_conf,iend_conf
203 call ene_eval(iprot,iii,ii)
205 call dawrite_ene(iprot,istart_conf,iend_conf,ientout)
209 do iii=istart_conf,iend_conf
210 write (icsa_rbank,vormat,rec=iii)
211 & iii,list_conf(iii,iprot),
212 & e_total(iii,iprot),
213 & ((nu(l,k,iii-istart_conf+1,iprot),l=1,natdim(k,iprot)),
214 & k=1,natlike(iprot))
220 if (.not.init_ene .and. mod_fourier(nloctyp).eq.0
221 & .and. .not. mod_elec .and. .not. mod_scp
222 & .and. .not. mod_tor .and. .not. mod_sccor) then
224 if (natlike(iprot).gt.0) close(icbase)
228 if (init_ene .or. mod_fourier(nloctyp).gt.0 .or. mod_elec
229 & .or. mod_scp .or. mod_tor .or. mod_sccor) close(ientout)
233 c print *,"Processor",me,me1," finished computing energies"
235 do i=1,scount(me1,iprot)
236 e_total_(i)=e_total(indstart(me1,iprot)+i-1,iprot)
237 eini_(i)=eini(indstart(me1,iprot)+i-1,iprot)
238 entfac_(i)=entfac(indstart(me1,iprot)+i-1,iprot)
239 rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot)
241 if (What.eq.2 .or. What.eq.3) then
242 c call MPI_Gatherv(e_total(indstart(me1,iprot),iprot),
243 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
244 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
246 call MPI_Gatherv(e_total_(1),
247 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
248 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
250 c call MPI_Gatherv(entfac(indstart(me1,iprot),iprot),
251 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot),
252 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
254 call MPI_Gatherv(entfac_(1),
255 & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot),
256 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
258 c call MPI_Gatherv(eini(indstart(me1,iprot),iprot),
259 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot),
260 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
262 call MPI_Gatherv(eini_(1),
263 & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot),
264 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
266 c call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
267 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
268 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
270 call MPI_Gatherv(rmstb_(1),
271 & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
272 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
276 c call MPI_Allgatherv(e_total(indstart(me1,iprot),iprot),
277 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
278 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
280 call MPI_Allgatherv(e_total_(1),
281 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
282 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
287 c Mean and lowest energies of structural classes
289 IF (NCHUNK_CONF(IPROT).EQ.1) THEN
292 do i=1,ntot_work(iprot)
296 do i=indstart(me1,iprot),indend(me1,iprot)
300 istart_conf=indstart(me1,iprot)
301 iend_conf=indend(me1,iprot)
303 do i=1,ntot_work(iprot)
307 iend_conf=ntot_work(iprot)
310 do i=1,ntot_work(iprot)
311 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
315 open (ientin,file=benefiles(iprot),status="old",
316 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
317 c call daread_ene(iprot,istart_conf,iend_conf)
318 call emin_search(iprot)
322 open (ientin,file=benefiles(iprot),status="old",
323 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
326 do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
327 iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
329 do istart_conf=1,ntot_work(iprot),maxstr_proc
330 iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
333 c Read the chunk of energies and derivatives off a DA scratchfile.
335 ipass_conf=ipass_conf+1
336 do i=1,ntot_work(iprot)
340 do i=istart_conf,iend_conf
345 write (iout,*) "ipass_conf",ipass_conf,
346 & " istart_conf",istart_conf," iend_conf",iend_conf
347 do i=1,ntot_work(iprot)
348 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
352 call daread_ene(iprot,istart_conf,iend_conf)
353 call emin_search(iprot)
362 c Complete the calculation of the lowest energies over all classes and
363 c distribute the values to all procs
366 write (iout,*) "Processor",me,me1," calling MPI_Reduce: 2"
367 do k=1,natlike(iprot)+2
368 write (iout,*) "iprot",iprot," k",k
369 do ib=1,nbeta(k,iprot)
370 write (iout,*) "ib",ib
371 write (iout,'(7hELOWEST,10e15.3)')
372 & elowest(ib,k,iprot)
377 do ibatch=1,natlike(iprot)+2
378 do ib=1,nbeta(ibatch,iprot)
379 elowest_aux(1,ib)=elowest(ib,ibatch,iprot)
380 elowest_aux(2,ib)=ind_lowest(ib,ibatch,iprot)
382 call MPI_Allreduce(elowest_aux(1,1),elowest_t(1,1),
383 & nbeta(ibatch,iprot),
384 & MPI_2DOUBLE_PRECISION, MPI_MINLOC, Comm1, IERROR)
386 do ib=1,nbeta(ibatch,iprot)
387 write (iout,*) "beta=",betaT(ib,ibatch,iprot)
388 write (iout,'(9helowest_t,10f15.3)')
389 & elowest_t(1,ib),elowest_t(2,ib)
391 write (iout,*) "Processor",me,me1," finished MPI_Reduce: 2"
393 do ib=1,nbeta(ibatch,iprot)
394 c write (iout,*) "IB",ib," ibatch",ibatch," iprot",iprot
396 elowest(ib,ibatch,iprot)=elowest_t(1,ib)
397 c write (iout,*) "elowest",ib,ibatch,iprot,
398 c & elowest(ib,ibatch,iprot)
400 ind_lowest(ib,ibatch,iprot)=elowest_t(2,ib)
408 c write (iout,*) "After averages, calling Thermal"
410 if (What.eq.2) call thermal(iprot)
411 c write (iout,*) "After thermal"
415 if (me1.eq.Master) then
419 write (iout,*) "iprot",iprot," elowest",
420 & ((elowest(ib,k,iprot),ib=1,nbeta(k,iprot)),k=1,natlike(iprot)+2)
423 do ib=1,nbeta(2,iprot)
424 fac=betaT(ib,2,iprot)
425 c 6/15/05 AL Heat capacity
426 zvtot(ib,iprot)=(esquare_ftot(ib,iprot)-
427 & emean_ftot(ib,iprot)**2)*fac*fac*Rgas
428 & -ebis_ftot(ib,iprot)+heatbase(iprot)
430 write (iout,*) "iprot",iprot," ib",ib,
431 & " emeantot",emean_ftot(ib,iprot),
432 & " esquaretot",esquare_ftot(ib,iprot),
433 & " ebistot",ebis_ftot(ib,iprot),
434 & " zvtot",zvtot(ib,iprot)
437 write (iout,*) "beta",
439 & " efree",efree(ib,2,iprot),
440 & " emean",emean_ftot(ib,iprot),
441 & " variance",esquare_ftot(ib,iprot)
444 c 4/13/04 AL Correlation coefficients between energy and nativelikeness
445 c in the respective classes.
447 write (iout,*) "iprot",iprot," ib",ib,
449 write (iout,*) "emean",emean_ftot(ib,iprot),
450 & " esquare",esquare_ftot(ib,iprot)
456 & "Averages of nativelikeness and energy-nativelikeness",
457 & " correlation coefficients:"
459 do i=1,natlike(iprot)
460 do ib=1,nbeta(i+1,iprot)
461 write (iout,*) "protein",iprot," property",i,
462 & " temperature",ib," natdim",natdim(i,iprot)," value(s)",
463 & (nuave(j,ib,i,iprot),j=1,natdim(i,iprot))
474 write (iout,*) "Calling CUTOFF_VIOLATION"
477 c write (iout,*) "CUTOFFEVAL",cutoffeval
479 if (cutoffeval) call cutoff_violation
480 t_func = t_func + tcpu() - tcpu_ini
482 write (iout,*) "Exitting FUNC1"
487 c-------------------------------------------------------------------------------
488 subroutine targetfun(n,x,nf,f,uiparm,ufparm)
491 include "DIMENSIONS.ZSCOPT"
494 integer IERROR,Errcode,status(MPI_STATUS_SIZE)
498 double precision x(n)
499 include "COMMON.WEIGHTS"
500 include "COMMON.XBOUND"
501 include "COMMON.IOUNITS"
502 include "COMMON.DERIV"
503 include "COMMON.VMCPAR"
504 include "COMMON.CLASSES"
506 include "COMMON.ENERGIES"
507 include "COMMON.WEIGHTDER"
508 include "COMMON.OPTIM"
509 include "COMMON.TIME1"
510 include "COMMON.COMPAR"
513 integer i,ii,it,inat,j,k,l,kk,iprot,ib
516 double precision urparm
518 double precision f,kara,viol,gnmr,gnmr1,aux
522 write (iout,*) me,me1,"Calling TARGETFUN"
524 call write_restart(n,x)
530 c Maximum likelihood contribution
532 do iT=1,nbeta(1,iprot)
533 f = f + weilik(iT,iprot)*sumlik(iT,iprot)
537 write (iout,*)"target_fun: sumlik"
539 write (iout,*) (sumlik(ii,iprot),ii=1,nbeta_l(iprot))
541 write (iout,*) "f before Cv =",f
545 c 6/15/05 AL Contribution from Cv
546 do ib=1,nbeta(2,iprot)
547 f=f+weiCv(ib,iprot)*(zvtot(ib,iprot)-target_cv(ib,iprot))**2
549 c write (iout,*) "target_fun from CV: f=",f
551 c write (iout,*) "target_fun: f=",f
552 c 8/2/2013 Contribution from conformation-dependent averages
554 do inat=1,natlike(iprot)
555 do ib=1,nbeta(inat+2,iprot)
556 do k=1,natdim(inat,iprot)
557 f = f+weinat(k,ib,inat,iprot)*
558 & (nuave(k,ib,inat,iprot)-nuexp(k,ib,inat,iprot))**2
560 write (iout,*) "iprot",iprot," inat",inat,
561 & " ib",ib," k=",k," nuave",nuave(k,ib,inat,iproti),
562 & " nuexp",nuexp(k,ib,inat,iprot),
563 & " weight",weinat(k,ib,inat,iprot)
569 c write (iout,*) "target_fun after adding nativelikeness: f=",f
571 c add constraints on weights
574 kara=kara+gnmr1(x(i),x_low(i),x_up(i))
579 c write (iout,*) "target_fun: f=",f,1.0d16*kara
583 c-----------------------------------------------------------------------------
584 subroutine targetgrad(n,x,nf,g,uiparm,rdum,ufparm)
587 include "DIMENSIONS.ZSCOPT"
588 include "COMMON.TIME1"
591 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
594 double precision x(n),
595 & g(n),gg(max_paropt)
596 include "COMMON.WEIGHTS"
597 include "COMMON.XBOUND"
598 include "COMMON.ENERGIES"
599 include "COMMON.CLASSES"
600 include "COMMON.IOUNITS"
601 include "COMMON.DERIV"
602 include "COMMON.COMPAR"
603 include "COMMON.WEIGHTDER"
604 include "COMMON.OPTIM"
605 include "COMMON.VMCPAR"
607 include "COMMON.NAMES"
614 double precision xi,delta /1.0d-6/
616 double precision urparm
618 double precision zewn(MaxT,maxprot),aux,zscder
619 double precision zewnvtot(maxT,maxprot)
620 double precision sumlik0(maxT,maxprot),zv0tot(maxT,maxprot),
621 & nuave0(maxdimnat,maxT,maxnatlike,maxprot)
622 double precision geps(nntyp),gtor(maxtor_var)
623 double precision rdum
625 integer ll1,ll2,nl1,nl2,nn
626 integer i,ii,ind,inat,j,k,kk,iprot,ib,n_weight
627 integer iti,itj,jj,icant
630 double precision gnmrprim,gnmr1prim
631 double precision tcpu_ini,tcpu
632 logical in_sc_pair_list
633 external in_sc_pair_list
639 write (iout,*) "Processor",me,me1," calling TARGETGRAD"
644 write (iout,*) "Calling TARGETGRAD"
649 write (iout,*) "x",(x(i),i=1,n)
656 do ib=1,nbeta(2,iprot)
657 c 6/15/05 AL Contributions from Cv
658 zewnvtot(ib,iprot)=2*weiCv(ib,iprot)*
659 & (zvtot(ib,iprot)-target_cv(ib,iprot))
671 if (mod_tor .or. mod_sccor) then
676 c Maximum likelihood gradient
678 do ib=1,nbeta(1,iprot)
681 if (imask(k).gt.0) then
683 g(kk)=g(kk)+weilik(ib,iprot)*sumlikder(kk,ib,iprot)
687 geps(k)=geps(k)+weilik(ib,iprot)*sumlikeps(k,ib,iprot)
690 gtor(k)=gtor(k)+weilik(ib,iprot)*sumliktor(k,ib,iprot)
695 write (iout,*) "gtor"
697 write(iout,*) k,gtor(k)
701 c Heat-capacity gradient
703 do ib=1,nbeta(2,iprot)
706 if (imask(k).gt.0) then
708 g(kk)=g(kk)+zvdertot(kk,ib,iprot)*zewnvtot(ib,iprot)
712 geps(k)=geps(k)+zewnvtot(ib,iprot)*zvepstot(k,ib,iprot)
715 gtor(k)=gtor(k)+zewnvtot(ib,iprot)*zvtortot(k,ib,iprot)
720 write (iout,*) "gtor"
722 write(iout,*) k,gtor(k)
724 write (iout,*) "grad: g",(g(i),i=1,n)
727 c Derivatives of conformational-dependent properties
729 do inat=1,natlike(iprot)
730 do ib=1,nbeta(inat+2,iprot)
731 call propderiv(ib,inat,iprot)
732 do l=1,natdim(inat,iprot)
733 aux=2*weinat(l,ib,inat,iprot)*
734 & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
736 write (iout,*) "iprot",iprot," ib",ib,
737 & " inat",inat," l",l,
738 & " nuave",nuave(l,ib,inat,iprot),
739 & " weight",weinat(ib,inat,iprot),
745 if (imask(k).gt.0) then
747 g(kk)=g(kk)+aux*rder(kk,l)
751 geps(k)=geps(k)+aux*reps(k,l)
754 gtor(k)=gtor(k)+aux*rtor(k,l)
761 write (iout,*) "gtor"
763 write(iout,*) k,gtor(k)
766 c Compute numerically the remaining part of the gradient
769 if (imask(i).gt.0) n_weight=n_weight+1
772 c Sum up the gradient in SC epsilons, if needed
777 g(ii)=g(ii)+geps(icant(iti,iti))
779 if (.not. in_sc_pair_list(iti,j)) then
780 g(ii)=g(ii)+0.5d0*geps(icant(iti,j))
788 g(ii)=geps(icant(iti,itj))
791 c write (iout,*) "n_ene",n_ene
792 if (mod_tor .or. mod_sccor) then
800 write (iout,*) "n",n," n_weight",n_weight
801 write (iout,*) "After SC: grad:"
803 write (iout,*) i,g(i)
809 do ib=1,nbeta(1,iprot)
810 sumlik0(ib,iprot)=sumlik(ib,iprot)
814 do ib=1,nbeta(2,iprot)
815 zv0tot(ib,iprot)=zvtot(ib,iprot)
819 do i=1,natlike(iprot)
820 do ib=1,nbeta(i+2,iprot)
821 do k=1,natdim(i,iprot)
822 nuave0(k,ib,i,iprot)=nuave(k,ib,i,iprot)
829 if (nbeta(2,iprot).gt.0) nn=nn-1
833 if (nbeta(2,iprot).gt.0) then
835 do ib=1,nbeta(2,iprot)
836 g(ii)=g(ii)+zewnvtot(ib,iprot)
841 write (iout,*) "n",n," n_weight",n_weight
842 write (iout,*) "After heat grad:"
844 write (iout,*) i,g(i)
848 c write (iout,*) "n_weight",n_weight," n",n
852 c write (iout,*) "i",i," x",x(i),xi
857 do j=1,nbeta(1,iprot)
858 g(i)=g(i)+weilik(ib,iprot)*(sumlik(j,iprot)
859 & -sumlik0(j,iprot))/delta
863 write (iout,*) "init numerical derivatives"
866 c Contribution from Cv
868 do ib=1,nbeta(2,iprot)
869 g(i)=g(i)+zewnvtot(ib,iprot)*
870 & (zvtot(ib,iprot)-zv0tot(ib,iprot))/delta
873 c Contribution from average nativelikenss and correlation coefficients
874 c between energy and nativelikeness
876 do ii=1,natlike(iprot)
877 do ib=1,nbeta(ii+2,iprot)
878 do k=1,natdim(i,iprot)
879 aux=2*weinat(k,ib,inat,iprot)*
880 & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
881 g(i)=g(i)+aux*(nuave(k,ib,ii,iprot)-
882 & nuave0(k,ib,ii,iprot))/delta
890 write (iout,*) "n",n," n_weight",n_weight
891 write (iout,*) "After num grad:"
893 write (iout,*) k,g(k)
899 write (iout,*) "finish numerical derivatives"
903 c 1/10/2007 AL Correction: ENDIF...IF insterted to prevent unmatched
904 c SEND is I_LOCAL_CHECK.GT.0
907 write (iout,*) "mode",mode
910 c transform the gradient to the x space
912 g(i)=g(i)*(x_up(i)-x_low(i))/(pi*(1.0d0+(x(i)
913 & -0.5d0*(x_up(i)+x_low(i)))**2))
916 c add constraints on weights
917 c write (iout,*) "Gradient in constraints"
919 g(i)=g(i)+1.0d16*gnmr1prim(x(i),x_low(i),x_up(i))
921 write (iout,*) i,x(i),x_low(i),x_up(i),
922 & gnmr1prim(x(i),x_low(i),x_up(i))
927 write (iout,*) "grad"
929 write (iout,*) i,g(i)
932 t_grad = t_grad + tcpu() - tcpu_ini
934 write (iout,*) "Exit targetgrad"
939 c------------------------------------------------------------------------------
940 subroutine matout(ndim,n,a)
941 implicit real*8 (a-h,o-z)
942 include "DIMENSIONS.ZSCOPT"
943 include "COMMON.IOUNITS"
944 double precision a(ndim)
945 character*6 line6 / '------' /
946 character*12 line12 / '------------' /
950 write (iout,1000) (k,k=i,nlim)
951 write (iout,1020) line6,(line12,k=i,nlim)
952 1000 format (/7x,6(i7,5x))
953 1020 format (a6,6a12)
958 3 b(k-i+1)=a(icant(i,j))
959 write (iout,1010) j,(b(k),k=1,jlim1)
962 1010 format (i3,3x,6(1pe11.4,1x))
965 c------------------------------------------------------------------------------
969 include "DIMENSIONS.ZSCOPT"
972 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
975 include "COMMON.WEIGHTS"
976 include "COMMON.ENERGIES"
977 include "COMMON.IOUNITS"
978 include "COMMON.DERIV"
979 include "COMMON.WEIGHTDER"
980 include "COMMON.VMCPAR"
982 include "COMMON.NAMES"
984 include "COMMON.CLASSES"
985 include "COMMON.THERMAL"
989 include "COMMON.TIME1"
990 double precision aux,zscder
991 double precision fac,fac2
992 double precision efj(2),emj(2),varj(2),edecum,edecum1,ebisdecum,
994 double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
995 & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
996 & enepsmixsqfj(nntyp,2),eavpfjprim(max_ene,2),
997 & entoravefj(0:3,maxtor_var,2),entoravefjprim(0:3,maxtor_var,2),
998 & entormixfj(0:3,maxtor_var,2),entormixsqfj(0:3,maxtor_var,2)
999 integer i,ii,j,k,kk,iprot,ib
1000 integer iti,itj,jj,icant
1006 do ib=1,nbeta(2,iprot)
1007 fac = betaT(ib,2,iprot)
1011 if (imask(k).gt.0) then
1013 ebisdecum=emix_pfbistot(k,ib,iprot)-
1014 & eave_pftot(k,ib,iprot)
1015 & *ebis_ftot(ib,iprot)
1016 edecum=emix_pfprimtot(k,ib,iprot)
1017 & -eave_pfprimtot(k,ib,iprot)
1018 & *emean_ftot(ib,iprot)
1019 edecum1=emix_pftot(k,ib,iprot)
1020 & -eave_pftot(k,ib,iprot)
1021 & *emean_ftot(ib,iprot)
1022 e2decum=emixsq_pftot(k,ib,iprot)
1023 & -eave_pftot(k,ib,iprot)
1024 & *esquare_ftot(ib,iprot)
1025 zvdertot(kk,ib,iprot)=Rgas*fac2*(2*edecum
1026 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1028 & -eave_pfbistot(k,ib,iprot)
1031 write (iout,*) "k=",k,
1032 & "eave_pftot",eave_pftot(k,ib,iprot),
1033 & " emix_pftot",emix_pftot(k,ib,iprot),
1034 & " emix_pfprimtot",emix_pfprimtot(k,ib,iprot),
1035 & " emix_pfbistot",emix_pfbistot(k,ib,iprot),
1036 & " eave_pfprimtot",eave_pfprimtot(k,ib,iprot),
1037 & " eave_pfbistot",eave_pfbistot(k,ib,iprot),
1038 & " emixsq_pftot",emixsq_pftot(k,ib,iprot),
1039 & " ebis_ftot",ebis_ftot(ib,iprot),
1040 & " emean_ftot",emean_ftot(ib,iprot)
1041 write (iout,*) "k=",k," edecum",edecum,
1042 & " edecum1",edecum1,
1043 & " e2decum",e2decum," ebisdecum",ebisdecum,
1045 & zvdertot(kk,ib,iprot)
1051 ebisdecum=eneps_mix_fbistot(k,ib,iprot)-
1052 & enepsave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1053 edecum=eneps_mix_ftot(k,ib,iprot)
1054 & -enepsave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1055 e2decum=eneps_mixsq_ftot(k,ib,iprot)
1056 & -enepsave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1057 zvepstot(k,ib,iprot)=ww(1)*(Rgas*fac2*(2*edecum
1058 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1059 & *edecum))+fac*ebisdecum)
1061 write (iout,*) ib,k,"ebisdecum",ebisdecum,
1062 & " edecum",edecum," e2decum",e2decum," zvepstot",
1063 & zvepstot(k,ibb,iprot)
1065 & " enepsave_ftot",enepsave_ftot(k,ib,iprot),
1066 & " eneps_mix_ftot",eneps_mix_ftot(k,ib,iprot),
1067 & " eneps_mix_fbistot",eneps_mix_fbistot(k,ib,iprot),
1068 & " eneps_mixsq_ftot",eneps_mixsq_ftot(k,ib,iprot),
1069 & " emean_ftot",emean_ftot(ib,iprot),
1070 & " ebis_ftot",ebis_ftot(ib,iprot)
1074 if (mod_tor .or. mod_sccor) then
1076 ebisdecum=entor_mix_fbistot(k,ib,iprot)-
1077 & entorave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1078 edecum=entor_mix_fprimtot(k,ib,iprot)
1079 & -entorave_fprimtot(k,ib,iprot)*emean_ftot(ib,iprot)
1080 edecum1=entor_mix_ftot(k,ib,iprot)
1081 & -entorave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1082 e2decum=entor_mixsq_ftot(k,ib,iprot)
1083 & -entorave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1084 zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1085 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1086 & *edecum1))-entorave_fbistot(k,ib,iprot)
1090 & " entorave_ftot",entorave_ftot(k,ib,iprot),
1091 & " entor_mix_ftot",entor_mix_ftot(k,ib,iprot),
1092 & " entor_mix_fprimtot",entor_mix_fprimtot(k,ib,iprot),
1093 & " entor_mix_fbistot",entor_mix_fbistot(k,ib,iprot),
1094 & " entorave_fprimtot",entorave_fprimtot(k,ib,iprot),
1095 & " entorave_fbistot",entorave_fbistot(k,ib,iprot),
1096 & " entor_mixsq_ftot",entor_mixsq_ftot(k,ib,iprot),
1097 & " emean_ftot",emean_ftot(ib,iprot),
1098 & " ebis_ftot",ebis_ftot(ib,iprot)
1099 write (iout,*) ib,k," ebisdecum",ebisdecum,
1100 & " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1101 & " zvtortot",zvtortot(k,ibb,iprot)
1106 write (iout,*) "iprot",iprot,
1108 write (iout,*) "zvder",
1109 & (zvdertot(k,ibb,iprot),k=1,kk)
1115 c------------------------------------------------------------------------------
1116 subroutine propderiv(ib,inat,iprot)
1117 c Compute the derivatives of conformation-dependent averages in force-field
1120 include "DIMENSIONS"
1121 include "DIMENSIONS.ZSCOPT"
1124 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1127 include "COMMON.WEIGHTS"
1128 include "COMMON.ENERGIES"
1129 include "COMMON.CLASSES"
1130 include "COMMON.IOUNITS"
1131 include "COMMON.DERIV"
1132 include "COMMON.WEIGHTDER"
1133 include "COMMON.VMCPAR"
1134 include "COMMON.GEO"
1135 include "COMMON.OPTIM"
1136 include "COMMON.NAMES"
1137 include "COMMON.COMPAR"
1139 include "COMMON.MPI"
1141 include "COMMON.TIME1"
1142 include "COMMON.VAR"
1143 double precision aux,zscder
1144 double precision fac
1145 double precision efj(2),emj(2),varj(2)
1146 double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1147 & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1148 & enepsmixsqfj(nntyp,2)
1150 integer i,ii,ind,j,k,kk,iprot,ib,ibb,l1
1151 integer iti,itj,jj,icant
1154 c write (iout,*) "Entered PROPDERIV"
1156 fac = betaT(ib,inat+2,iprot)
1157 c write (iout,*) "derivatives: iprot",iprot,
1159 do l=1,natdim(inat,iprot)
1162 if (imask(k).gt.0) then
1164 rder(kk,l)=-fac*(nu_pf(k,l,ib,inat,iprot)
1165 & -nuave(l,ib,inat,iprot)*
1166 & eave_nat_pftot(k,ib,inat,iprot))
1171 do l=1,natdim(inat,iprot)
1173 reps(k,l)=-ww(1)*fac*(nuepsave_f(k,l,ib,inat,iprot)-
1174 & nuave(l,ib,inat,iprot)*enepsave_nat_ftot(k,ib,inat,iprot))
1178 if (mod_tor .or. mod_sccor) then
1179 do l=1,natdim(inat,iprot)
1181 rtor(k,l)=-fac*(nutorave_f(k,l,ib,inat,iprot)-
1182 & nuave(l,inat,ib,iprot)*entorave_nat_ftot(k,ib,inat,iprot))
1188 c-------------------------------------------------------------------------------
1189 subroutine ene_eval(iprot,i,ii)
1191 include "DIMENSIONS"
1192 include "DIMENSIONS.ZSCOPT"
1195 integer IERROR,ErrCode
1196 include "COMMON.MPI"
1198 include "COMMON.CONTROL"
1199 include "COMMON.COMPAR"
1200 include "COMMON.WEIGHTS"
1201 include "COMMON.WEIGHTDER"
1202 include "COMMON.CLASSES"
1203 include "COMMON.ENERGIES"
1204 include "COMMON.IOUNITS"
1205 include "COMMON.VMCPAR"
1206 include "COMMON.NAMES"
1207 include "COMMON.INTERACT"
1208 include "COMMON.TIME1"
1209 include "COMMON.CHAIN"
1210 include "COMMON.VAR"
1211 include "COMMON.GEO"
1212 include "COMMON.SCCOR"
1213 include "COMMON.TORSION"
1214 C Define local variables
1215 integer i,ii,ij,jj,j,k,ik,jk,l,ll,lll,iprot,iblock
1216 double precision energia(0:max_ene)
1217 double precision etoti,enepsjk
1218 double precision ftune_eps
1226 write (iout,*) "ene_eval: iprot",iprot," i",i," ii",ii
1228 if (init_ene .or. (mod_fourier(nloctyp).gt.0
1229 & .or. mod_elec .or. mod_scp .or. mod_tor .or. mod_sccor) ) then
1231 write (iout,*) "FUNC1: doing UNRES"
1232 write (iout,*) "Protein",iprot," i",i," ii",ii," nres",nres
1235 call restore_ccoords(iprot,ii)
1236 call int_from_cart1(.false.)
1238 write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1239 write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1240 write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1241 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1242 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1243 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1244 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1245 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1250 call etotal(energia(0))
1251 c if (natlike(iprot).gt.0)
1252 c & call natsimil(i,ii,iprot,.false.)
1253 e_total(i,iprot)=energia(0)
1255 enetb(ii,j,iprot)=energia(j)
1258 write (iout,*) "FUN nntyp",nntyp," eneps",eneps_temp(1,45)
1261 eneps(k,j,ii,iprot)=eneps_temp(k,j)
1265 write (iout,'(2i5,18(1pe12.4))') i,ii,
1266 & (enetb(ii,j,iprot),j=1,n_ene)
1267 write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1268 & i,energia(0),eini(i,iprot),entfac(i,iprot)
1271 c write (iout,*) "Matrix eneps"
1274 c write (iout,*) j,k,icant(j,k),
1275 c & eneps_temp(1,icant(j,k)),
1276 c & eneps_temp(2,icant(j,k))
1282 if (mod_tor .or. mod_sccor) then
1283 c write (iout,*) "etor",enetb(ii,13,iprot),enetb(ii,19,iprot)
1286 do ij=-ntortyp+1,ntortyp-1
1287 do jj=-ntortyp+1,ntortyp-1
1288 if (mask_tor(0,jj,ij,iblock).eq.1) then
1290 entor(k,ii,iprot)=etor_temp(0,0,jj,ij,iblock)
1291 c write (iout,*) "etor_temp",restyp(itortyp(j)),
1292 c & restyp(itortyp(i)),
1293 c & k,etor_temp(ji,ii)
1294 else if (mask_tor(0,jj,ij,iblock).gt.1) then
1295 if (tor_mode.eq.0) then
1296 do l=1,nterm(jj,ij,iblock)
1298 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1300 if (mask_tor(0,jj,ij,iblock).gt.2) then
1301 do l=nterm(jj,ij,iblock)+1,2*nterm(jj,ij,iblock)
1303 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1306 else if (tor_mode.eq.1) then
1308 do l=1,nterm_kcc(jj,ij)
1309 do ll=1,nterm_kcc_Tb(jj,ij)
1310 do lll=1,nterm_kcc_Tb(jj,ij)
1312 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1316 if (mask_tor(0,jj,ij,1).gt.2) then
1317 do l=nterm_kcc(jj,ij)+1,2*nterm_kcc(jj,ij)
1318 do ll=1,nterm_kcc_Tb(jj,ij)
1319 do lll=1,nterm_kcc_Tb(jj,ij)
1321 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1332 do ij=-nsccortyp+1,nsccortyp-1
1333 do jj=-nsccortyp+1,nsccortyp-1
1335 if (mask_tor(ll,jj,ij,1).eq.1) then
1337 entor(k,ii,iprot)=etor_temp(0,ll,jj,ij,1)
1338 c write (iout,*) "etor_temp",restyp(isccortyp(jj)),
1339 c & restyp(isccortyp(ij)),
1340 c & k,etor_temp(jj,ij)
1341 else if (mask_tor(ll,jj,ij,1).gt.1) then
1342 do l=1,nterm_sccor(jj,ij)
1344 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1346 do l=nterm_sccor(jj,ij)+1,2*nterm_sccor(jj,ij)
1348 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1356 c write (iout,*) "Conformation",i
1357 c write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1358 c write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1359 c write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1360 c write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1361 c write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1362 c write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1363 c write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1364 c write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1365 c call enerprint(energia(0))
1366 write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1367 & i,energia(0),eini(i,iprot),entfac(i,iprot)
1368 write (iout,'(i5,18(1pe12.4))') ii,
1369 & (enetb(ii,j,iprot),j=1,n_ene)
1370 c write (iout,'(i5,16(1pe12.4))') i,
1371 c & (energia(j),j=1,n_ene),energia(0)
1372 c print *,"Processor",me,me1," computing energies",i
1374 if (energia(0).ge.1.0d99) then
1375 if (me1.eq.Master) then
1376 write (iout,*) "FUNC1:CHUJ NASTAPIL in energy evaluation for",
1377 & " point",i,". Probably NaNs in some of the energy components."
1378 write (iout,*) "The components of the energy are:"
1379 call enerprint(energia(0))
1380 write (iout,*) "Conformation",i,ii
1381 write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1382 write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1383 write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1384 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1385 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1386 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1387 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1388 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1389 write (iout,*) "Calculation terminated at this point.",
1390 & " Check the database of conformations"
1394 call MPI_Abort(ALL_COMM,ErrCode,IERROR)
1396 stop "SEVERE error in energy calculation"
1400 enetb(ii,1,iprot)=0.0d0
1403 enetb(ii,1,iprot)=enetb(ii,1,iprot)+
1404 & ftune_eps(eps(j,k))*eneps(1,icant(j,k),ii,iprot)+
1405 & eps(j,k)*eneps(2,icant(j,k),ii,iprot)
1411 etoti=etoti+enetb(ii,j,iprot)*ww(j)
1413 e_total(i,iprot)=etoti
1414 c write (iout,*) "FUNC1 natlike",iprot,natlike(iprot)
1415 if (natlike(iprot).gt.0) then
1416 call restore_ccoords(iprot,ii)
1417 call int_from_cart1(.false.)
1418 c call natsimil(i,ii,iprot,.false.)
1420 c & "natsimil",i,ii,(nu(k,ii,iprot),k=1,natlike(iprot))
1423 write (iout,'(i5,18(1pe12.4))') i,
1424 & (enetb(ii,j,iprot),j=1,n_ene)
1425 write (iout,'(18(1pe12.4))')
1426 & (eneps(1,j,ii,iprot),j=201,210)
1427 write(iout,'(4hENER,i10,3f15.3)')
1428 & i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)