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
58 write (iout,*) "FUNC1: Init_Ene ",init_ene
60 write (iout,*) "ELOWEST at the beginning of FUC1"
62 do ibatch=1,natlike(iprot)+2
63 do ib=1,nbeta(ibatch,iprot)
64 write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
65 & " elowest",elowest(ib,ibatch,iprot)
71 call MPI_Bcast( What, 1, MPI_INTEGER, Master, Comm1, IERROR)
72 call MPI_Bcast( cutoffeval, 1, MPI_LOGICAL, Master, Comm1, IERROR)
73 call MPI_Bcast(x(1),n,MPI_DOUBLE_PRECISION, Master, Comm1, IERROR)
76 write (iout,*) "Processor",me,me1," What",What
77 write (iout,*) "Processor",me,me1," x",(x(i),i=1,n)
82 write (iout,*) "x",(x(i),i=1,n)
86 C Convert variables into weights and other UNRES energy parameters.
88 C Calculate the total energies.
90 c write (iout,*) "Processor",me," nprot",nprot
91 c Revise the list of conformations
93 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 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," mod_ang",mod_ang
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 c Change AL 12/30/2017
161 if (.not. mod_other_params)
162 & open (ientout,file=benefiles(iprot),status="unknown",
163 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
167 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
169 iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
171 do i=1,ntot_work(iprot),maxstr_proc
173 iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
177 c Read the chunk of conformations off a DA scratchfile.
179 if (.not. init_ene .and. mod_fourier(nloctyp).eq.0
180 & .and. .not. mod_elec .and. .not. mod_scp
181 & .and. .not. mod_tor .and. .not. mod_sccor .and. .not. mod_ang)
184 write (iout,*) "func1: calling daread_ene",istart_conf,
187 call daread_ene(iprot,istart_conf,iend_conf)
188 if (natlike(iprot).gt.0) then
190 write (iout,*) "Calling daread_coords",istart_conf,
193 call daread_ccoords(iprot,istart_conf,iend_conf)
195 do iii=istart_conf,iend_conf
196 call ene_eval(iprot,iii,ii)
199 write (iout,*) "func1: calling dawrite_ene",istart_conf,
202 call dawrite_ene(iprot,istart_conf,iend_conf,ientin)
204 call daread_ccoords(iprot,istart_conf,iend_conf)
205 do iii=istart_conf,iend_conf
206 call ene_eval(iprot,iii,ii)
208 c AL 12/30/17 - writing energies is not needed when whole energy is evaluated
209 c call dawrite_ene(iprot,istart_conf,iend_conf,ientout)
213 do iii=istart_conf,iend_conf
214 write (icsa_rbank,vormat,rec=iii)
215 & iii,list_conf(iii,iprot),
216 & e_total(iii,iprot),
217 & ((nu(l,k,iii-istart_conf+1,iprot),l=1,natdim(k,iprot)),
218 & k=1,natlike(iprot))
224 if (.not.init_ene .and. mod_fourier(nloctyp).eq.0
225 & .and. .not. mod_elec .and. .not. mod_scp
226 & .and. .not. mod_tor .and. .not. mod_sccor .and. .not. mod_ang)
229 if (natlike(iprot).gt.0) close(icbase)
233 if (init_ene .or. mod_fourier(nloctyp).gt.0 .or. mod_elec
234 & .or. mod_scp .or. mod_tor .or. mod_sccor .or. mod_ang)
239 c print *,"Processor",me,me1," finished computing energies"
241 do i=1,scount(me1,iprot)
242 e_total_(i)=e_total(indstart(me1,iprot)+i-1,iprot)
243 eini_(i)=eini(indstart(me1,iprot)+i-1,iprot)
244 entfac_(i)=entfac(indstart(me1,iprot)+i-1,iprot)
245 rmstb_(i)=rmstb(indstart(me1,iprot)+i-1,iprot)
247 if (What.eq.2 .or. What.eq.3) then
248 c call MPI_Gatherv(e_total(indstart(me1,iprot),iprot),
249 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
250 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
252 call MPI_Gatherv(e_total_(1),
253 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
254 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
256 c call MPI_Gatherv(entfac(indstart(me1,iprot),iprot),
257 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot),
258 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
260 call MPI_Gatherv(entfac_(1),
261 & scount(me1,iprot),MPI_DOUBLE_PRECISION,entfac(1,iprot),
262 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
264 c call MPI_Gatherv(eini(indstart(me1,iprot),iprot),
265 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot),
266 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
268 call MPI_Gatherv(eini_(1),
269 & scount(me1,iprot),MPI_DOUBLE_PRECISION,eini(1,iprot),
270 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
272 c call MPI_Gatherv(rmstb(indstart(me1,iprot),iprot),
273 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
274 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
276 call MPI_Gatherv(rmstb_(1),
277 & scount(me1,iprot),MPI_DOUBLE_PRECISION,rmstb(1,iprot),
278 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,Master,
282 c call MPI_Allgatherv(e_total(indstart(me1,iprot),iprot),
283 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
284 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
286 call MPI_Allgatherv(e_total_(1),
287 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total(1,iprot),
288 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
293 c Mean and lowest energies of structural classes
295 IF (NCHUNK_CONF(IPROT).EQ.1) THEN
298 do i=1,ntot_work(iprot)
302 do i=indstart(me1,iprot),indend(me1,iprot)
306 istart_conf=indstart(me1,iprot)
307 iend_conf=indend(me1,iprot)
309 do i=1,ntot_work(iprot)
313 iend_conf=ntot_work(iprot)
316 do i=1,ntot_work(iprot)
317 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
321 C Change AL 12/30/2017
322 if (.not.mod_other_params)
323 &open (ientin,file=benefiles(iprot),status="old",
324 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
325 c call daread_ene(iprot,istart_conf,iend_conf)
326 call emin_search(iprot)
330 if (.not.mod_other_params)
331 &open (ientin,file=benefiles(iprot),status="old",
332 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
335 do istart_conf=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
336 iend_conf=min0(istart_conf+maxstr_proc-1,indend(me1,iprot))
338 do istart_conf=1,ntot_work(iprot),maxstr_proc
339 iend_conf=min0(istart_conf+maxstr_proc-1,ntot_work(iprot))
342 c Read the chunk of energies and derivatives off a DA scratchfile.
344 ipass_conf=ipass_conf+1
345 do i=1,ntot_work(iprot)
349 do i=istart_conf,iend_conf
354 write (iout,*) "ipass_conf",ipass_conf,
355 & " istart_conf",istart_conf," iend_conf",iend_conf
356 do i=1,ntot_work(iprot)
357 write (iout,*) "i",i," i2ii",i2ii(i,iprot)
361 call daread_ene(iprot,istart_conf,iend_conf)
362 call emin_search(iprot)
371 c Complete the calculation of the lowest energies over all classes and
372 c distribute the values to all procs
375 write (iout,*) "Processor",me,me1," calling MPI_Reduce: 2"
376 do k=1,natlike(iprot)+2
377 write (iout,*) "iprot",iprot," k",k
378 do ib=1,nbeta(k,iprot)
379 write (iout,*) "ib",ib
380 write (iout,'(7hELOWEST,10e15.3)')
381 & elowest(ib,k,iprot)
386 do ibatch=1,natlike(iprot)+2
387 do ib=1,nbeta(ibatch,iprot)
388 elowest_aux(1,ib)=elowest(ib,ibatch,iprot)
389 elowest_aux(2,ib)=ind_lowest(ib,ibatch,iprot)
391 call MPI_Allreduce(elowest_aux(1,1),elowest_t(1,1),
392 & nbeta(ibatch,iprot),
393 & MPI_2DOUBLE_PRECISION, MPI_MINLOC, Comm1, IERROR)
395 do ib=1,nbeta(ibatch,iprot)
396 write (iout,*) "beta=",betaT(ib,ibatch,iprot)
397 write (iout,'(9helowest_t,10f15.3)')
398 & elowest_t(1,ib),elowest_t(2,ib)
400 write (iout,*) "Processor",me,me1," finished MPI_Reduce: 2"
402 do ib=1,nbeta(ibatch,iprot)
403 c write (iout,*) "IB",ib," ibatch",ibatch," iprot",iprot
405 elowest(ib,ibatch,iprot)=elowest_t(1,ib)
406 c write (iout,*) "elowest",ib,ibatch,iprot,
407 c & elowest(ib,ibatch,iprot)
409 ind_lowest(ib,ibatch,iprot)=elowest_t(2,ib)
417 c write (iout,*) "After averages, calling Thermal"
419 if (What.eq.2) call thermal(iprot)
420 c write (iout,*) "After thermal"
424 if (me1.eq.Master) then
428 write (iout,*) "iprot",iprot," elowest",
429 & ((elowest(ib,k,iprot),ib=1,nbeta(k,iprot)),k=1,natlike(iprot)+2)
432 do ib=1,nbeta(2,iprot)
433 fac=betaT(ib,2,iprot)
434 c 6/15/05 AL Heat capacity
435 zvtot(ib,iprot)=(esquare_ftot(ib,iprot)-
436 & emean_ftot(ib,iprot)**2)*fac*fac*Rgas
437 & -ebis_ftot(ib,iprot)+heatbase(iprot)
439 write (iout,*) "iprot",iprot," ib",ib,
440 & " emeantot",emean_ftot(ib,iprot),
441 & " esquaretot",esquare_ftot(ib,iprot),
442 & " ebistot",ebis_ftot(ib,iprot),
443 & " zvtot",zvtot(ib,iprot)
446 write (iout,*) "beta",
448 & " efree",efree(ib,2,iprot),
449 & " emean",emean_ftot(ib,iprot),
450 & " variance",esquare_ftot(ib,iprot)
453 c 4/13/04 AL Correlation coefficients between energy and nativelikeness
454 c in the respective classes.
456 write (iout,*) "iprot",iprot," ib",ib,
458 write (iout,*) "emean",emean_ftot(ib,iprot),
459 & " esquare",esquare_ftot(ib,iprot)
465 & "Averages of nativelikeness and energy-nativelikeness",
466 & " correlation coefficients:"
468 do i=1,natlike(iprot)
469 do ib=1,nbeta(i+1,iprot)
470 write (iout,*) "protein",iprot," property",i,
471 & " temperature",ib," natdim",natdim(i,iprot)," value(s)",
472 & (nuave(j,ib,i,iprot),j=1,natdim(i,iprot))
483 write (iout,*) "Calling CUTOFF_VIOLATION"
486 c write (iout,*) "CUTOFFEVAL",cutoffeval
488 if (cutoffeval) call cutoff_violation
489 t_func = t_func + tcpu() - tcpu_ini
491 write (iout,*) "Exitting FUNC1"
496 c-------------------------------------------------------------------------------
497 subroutine targetfun(n,x,nf,f,uiparm,ufparm)
500 include "DIMENSIONS.ZSCOPT"
503 integer IERROR,Errcode,status(MPI_STATUS_SIZE)
507 double precision x(n)
508 include "COMMON.WEIGHTS"
509 include "COMMON.IOUNITS"
510 include "COMMON.DERIV"
511 include "COMMON.VMCPAR"
512 include "COMMON.CLASSES"
514 include "COMMON.ENERGIES"
515 include "COMMON.WEIGHTDER"
516 include "COMMON.OPTIM"
517 include "COMMON.TIME1"
518 include "COMMON.COMPAR"
521 integer i,ii,it,inat,j,k,l,kk,iprot,ib
524 double precision urparm
526 double precision f,kara,viol,gnmr,gnmr1,aux
531 write (iout,*) me,me1,"Calling TARGETFUN"
532 write (iout,*) "n",n," x",(x(i),i=1,n)
535 call write_restart(n,x)
541 c Maximum likelihood contribution
543 do iT=1,nbeta(1,iprot)
544 f = f + weilik(iT,iprot)*sumlik(iT,iprot)
548 write (iout,*)"target_fun: sumlik"
550 write (iout,*) (sumlik(ii,iprot),ii=1,nbeta(1,iprot))
552 write (iout,*) "f before Cv =",f
556 c 6/15/05 AL Contribution from Cv
557 do ib=1,nbeta(2,iprot)
558 f=f+weiCv(ib,iprot)*(zvtot(ib,iprot)-target_cv(ib,iprot))**2
560 c write (iout,*) "target_fun from CV: f=",f
562 c write (iout,*) "target_fun: f=",f
563 c 8/2/2013 Contribution from conformation-dependent averages
565 do inat=1,natlike(iprot)
566 do ib=1,nbeta(inat+2,iprot)
567 do k=1,natdim(inat,iprot)
568 f = f+weinat(k,ib,inat,iprot)*
569 & (nuave(k,ib,inat,iprot)-nuexp(k,ib,inat,iprot))**2
571 write (iout,*) "iprot",iprot," inat",inat,
572 & " ib",ib," k=",k," nuave",nuave(k,ib,inat,iprot),
573 & " nuexp",nuexp(k,ib,inat,iprot),
574 & " weight",weinat(k,ib,inat,iprot)
580 c write (iout,*) "target_fun after adding nativelikeness: f=",f
582 c add constraints on weights
585 kara=kara+gnmr1(x(i),x_low(i),x_up(i))
591 write (iout,*) "target_fun: f=",f,1.0d16*kara
597 c-----------------------------------------------------------------------------
598 subroutine targetgrad(n,x,nf,g,uiparm,rdum,ufparm)
601 include "DIMENSIONS.ZSCOPT"
602 include "COMMON.TIME1"
605 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
608 double precision x(n),
609 & g(n),gg(max_paropt)
610 include "COMMON.WEIGHTS"
611 include "COMMON.ENERGIES"
612 include "COMMON.CLASSES"
613 include "COMMON.IOUNITS"
614 include "COMMON.DERIV"
615 include "COMMON.COMPAR"
616 include "COMMON.WEIGHTDER"
617 include "COMMON.OPTIM"
618 include "COMMON.VMCPAR"
620 include "COMMON.NAMES"
627 double precision xi,delta /1.0d-6/
629 double precision urparm
631 double precision zewn(MaxT,maxprot),aux,zscder
632 double precision zewnvtot(maxT,maxprot)
633 double precision sumlik0(maxT,maxprot),zv0tot(maxT,maxprot),
634 & nuave0(maxdimnat,maxT,maxnatlike,maxprot)
635 double precision geps(nntyp),gtor(maxtor_var),gang(maxang_var)
636 double precision rdum
638 integer ll1,ll2,nl1,nl2,nn
639 integer i,ii,ind,inat,j,k,kk,iprot,ib,n_weight
640 integer iti,itj,jj,icant
643 double precision gnmrprim,gnmr1prim
644 double precision tcpu_ini,tcpu
645 logical in_sc_pair_list
646 external in_sc_pair_list
647 integer ind_shield /25/
653 write (iout,*) "Processor",me,me1," calling TARGETGRAD"
658 write (iout,*) "Calling TARGETGRAD"
663 write (iout,*) "x",(x(i),i=1,n)
670 do ib=1,nbeta(2,iprot)
671 c 6/15/05 AL Contributions from Cv
672 zewnvtot(ib,iprot)=2*weiCv(ib,iprot)*
673 & (zvtot(ib,iprot)-target_cv(ib,iprot))
685 if (mod_tor .or. mod_sccor) then
695 c Maximum likelihood gradient
697 do ib=1,nbeta(1,iprot)
700 if (imask(k).gt.0 .and. k.ne.ind_shield) then
702 g(kk)=g(kk)+weilik(ib,iprot)*sumlikder(kk,ib,iprot)
706 geps(k)=geps(k)+weilik(ib,iprot)*sumlikeps(k,ib,iprot)
709 gtor(k)=gtor(k)+weilik(ib,iprot)*sumliktor(k,ib,iprot)
712 gang(k)=gang(k)+weilik(ib,iprot)*sumlikang(k,ib,iprot)
717 write (iout,*) "gtor"
719 write(iout,*) k,gtor(k)
721 write (iout,*) "nang_var",nang_var," gang"
723 write(iout,*) k,gang(k)
727 c Heat-capacity gradient
729 do ib=1,nbeta(2,iprot)
732 if (imask(k).gt.0 .and. k.ne.ind_shield) then
734 g(kk)=g(kk)+zvdertot(kk,ib,iprot)*zewnvtot(ib,iprot)
738 geps(k)=geps(k)+zewnvtot(ib,iprot)*zvepstot(k,ib,iprot)
741 gtor(k)=gtor(k)+zewnvtot(ib,iprot)*zvtortot(k,ib,iprot)
744 gang(k)=gang(k)+zewnvtot(ib,iprot)*zvangtot(k,ib,iprot)
749 write (iout,*) "gtor"
751 write(iout,*) k,gtor(k)
753 write (iout,*) "gang"
755 write(iout,*) k,gang(k)
757 write (iout,*) "grad: g",(g(i),i=1,n)
760 c Derivatives of conformational-dependent properties
762 do inat=1,natlike(iprot)
763 do ib=1,nbeta(inat+2,iprot)
764 call propderiv(ib,inat,iprot)
765 do l=1,natdim(inat,iprot)
766 aux=2*weinat(l,ib,inat,iprot)*
767 & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
769 write (iout,*) "iprot",iprot," ib",ib,
770 & " inat",inat," l",l,
771 & " nuave",nuave(l,ib,inat,iprot),
772 & " weight",weinat(ib,inat,iprot),
778 if (imask(k).gt.0 .and. k.ne.ind_shield) then
780 g(kk)=g(kk)+aux*rder(kk,l)
784 geps(k)=geps(k)+aux*reps(k,l)
787 gtor(k)=gtor(k)+aux*rtor(k,l)
790 gang(k)=gang(k)+aux*rang(k,l)
797 write (iout,*) "gtor"
799 write(iout,*) k,gtor(k)
801 write (iout,*) "gang"
803 write(iout,*) k,gang(k)
806 c Compute numerically the remaining part of the gradient
809 if (imask(i).gt.0 .and. i.ne.ind_shield) n_weight=n_weight+1
812 c Sum up the gradient in SC epsilons, if needed
817 g(ii)=g(ii)+geps(icant(iti,iti))
819 if (.not. in_sc_pair_list(iti,j)) then
820 g(ii)=g(ii)+0.5d0*geps(icant(iti,j))
828 g(ii)=geps(icant(iti,itj))
831 c write (iout,*) "n_ene",n_ene
832 if (mod_tor .or. mod_sccor) then
846 write (iout,*) "n",n," n_weight",n_weight
847 write (iout,*) "After SC: grad:"
849 write (iout,*) i,g(i)
855 do ib=1,nbeta(1,iprot)
856 sumlik0(ib,iprot)=sumlik(ib,iprot)
860 do ib=1,nbeta(2,iprot)
861 zv0tot(ib,iprot)=zvtot(ib,iprot)
865 do i=1,natlike(iprot)
866 do ib=1,nbeta(i+2,iprot)
867 do k=1,natdim(i,iprot)
868 nuave0(k,ib,i,iprot)=nuave(k,ib,i,iprot)
875 if (nbeta(2,iprot).gt.0) nn=nn-1
879 if (nbeta(2,iprot).gt.0) then
881 do ib=1,nbeta(2,iprot)
882 g(ii)=g(ii)+zewnvtot(ib,iprot)
887 write (iout,*) "n",n," n_weight",n_weight
888 write (iout,*) "After heat grad:"
890 write (iout,*) i,g(i)
894 c write (iout,*) "n_weight",n_weight," n",n
898 c write (iout,*) "i",i," x",x(i),xi
903 do ib=1,nbeta(1,iprot)
904 g(i)=g(i)+weilik(ib,iprot)*(sumlik(ib,iprot)
905 & -sumlik0(ib,iprot))/delta
906 c write (iout,*) "iprot",iprot," g",g(i)
910 write (iout,*) "init numerical derivatives"
913 c Contribution from Cv
915 do ib=1,nbeta(2,iprot)
916 g(i)=g(i)+zewnvtot(ib,iprot)*
917 & (zvtot(ib,iprot)-zv0tot(ib,iprot))/delta
920 c Contribution from average nativelikenss and correlation coefficients
921 c between energy and nativelikeness
923 do ii=1,natlike(iprot)
924 do ib=1,nbeta(ii+2,iprot)
925 do k=1,natdim(i,iprot)
926 aux=2*weinat(k,ib,inat,iprot)*
927 & (nuave(l,ib,inat,iprot)-nuexp(l,ib,inat,iprot))
928 g(i)=g(i)+aux*(nuave(k,ib,ii,iprot)-
929 & nuave0(k,ib,ii,iprot))/delta
937 write (iout,*) "n",n," n_weight",n_weight
938 write (iout,*) "After num grad:"
940 write (iout,*) k,g(k)
946 write (iout,*) "finish numerical derivatives"
950 c 1/10/2007 AL Correction: ENDIF...IF insterted to prevent unmatched
951 c SEND is I_LOCAL_CHECK.GT.0
954 write (iout,*) "mode",mode
957 c transform the gradient to the x space
959 g(i)=g(i)*(x_up(i)-x_low(i))/(pi*(1.0d0+(x(i)
960 & -0.5d0*(x_up(i)+x_low(i)))**2))
963 c add constraints on weights
964 c write (iout,*) "Gradient in constraints"
966 g(i)=g(i)+1.0d16*gnmr1prim(x(i),x_low(i),x_up(i))
968 write (iout,*) i,x(i),x_low(i),x_up(i),
969 & gnmr1prim(x(i),x_low(i),x_up(i))
974 write (iout,*) "grad"
976 write (iout,*) i,g(i)
979 t_grad = t_grad + tcpu() - tcpu_ini
981 write (iout,*) "Exit targetgrad"
986 c------------------------------------------------------------------------------
987 subroutine matout(ndim,n,a)
988 implicit real*8 (a-h,o-z)
989 include "DIMENSIONS.ZSCOPT"
990 include "COMMON.IOUNITS"
991 double precision a(ndim)
992 character*6 line6 / '------' /
993 character*12 line12 / '------------' /
997 write (iout,1000) (k,k=i,nlim)
998 write (iout,1020) line6,(line12,k=i,nlim)
999 1000 format (/7x,6(i7,5x))
1000 1020 format (a6,6a12)
1005 3 b(k-i+1)=a(icant(i,j))
1006 write (iout,1010) j,(b(k),k=1,jlim1)
1009 1010 format (i3,3x,6(1pe11.4,1x))
1012 c------------------------------------------------------------------------------
1015 include "DIMENSIONS"
1016 include "DIMENSIONS.ZSCOPT"
1019 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1022 include "COMMON.WEIGHTS"
1023 include "COMMON.ENERGIES"
1024 include "COMMON.IOUNITS"
1025 include "COMMON.DERIV"
1026 include "COMMON.WEIGHTDER"
1027 include "COMMON.VMCPAR"
1028 include "COMMON.GEO"
1029 include "COMMON.NAMES"
1030 include "COMMON.VAR"
1031 include "COMMON.CLASSES"
1032 include "COMMON.THERMAL"
1034 include "COMMON.MPI"
1036 include "COMMON.TIME1"
1037 double precision aux,zscder
1038 double precision fac,fac2
1039 double precision efj(2),emj(2),varj(2),edecum,edecum1,ebisdecum,
1041 double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1042 & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1043 & enepsmixsqfj(nntyp,2),eavpfjprim(max_ene,2),
1044 & entoravefj(0:3,maxtor_var,2),entoravefjprim(0:3,maxtor_var,2),
1045 & entormixfj(0:3,maxtor_var,2),entormixsqfj(0:3,maxtor_var,2),
1046 & enangavefj(0:3,maxang_var,2),enangavefjprim(0:3,maxang_var,2),
1047 & enangmixfj(0:3,maxang_var,2),enangmixsqfj(0:3,maxang_var,2)
1048 integer i,ii,j,k,kk,iprot,ib
1049 integer iti,itj,jj,icant
1052 integer ind_shield /25/
1056 do ib=1,nbeta(2,iprot)
1057 fac = betaT(ib,2,iprot)
1061 if (imask(k).gt.0 .and. k.ne.ind_shield) then
1063 ebisdecum=emix_pfbistot(k,ib,iprot)-
1064 & eave_pftot(k,ib,iprot)
1065 & *ebis_ftot(ib,iprot)
1066 edecum=emix_pfprimtot(k,ib,iprot)
1067 & -eave_pfprimtot(k,ib,iprot)
1068 & *emean_ftot(ib,iprot)
1069 edecum1=emix_pftot(k,ib,iprot)
1070 & -eave_pftot(k,ib,iprot)
1071 & *emean_ftot(ib,iprot)
1072 e2decum=emixsq_pftot(k,ib,iprot)
1073 & -eave_pftot(k,ib,iprot)
1074 & *esquare_ftot(ib,iprot)
1075 zvdertot(kk,ib,iprot)=Rgas*fac2*(2*edecum
1076 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1078 & -eave_pfbistot(k,ib,iprot)
1081 write (iout,*) "k=",k,
1082 & "eave_pftot",eave_pftot(k,ib,iprot),
1083 & " emix_pftot",emix_pftot(k,ib,iprot),
1084 & " emix_pfprimtot",emix_pfprimtot(k,ib,iprot),
1085 & " emix_pfbistot",emix_pfbistot(k,ib,iprot),
1086 & " eave_pfprimtot",eave_pfprimtot(k,ib,iprot),
1087 & " eave_pfbistot",eave_pfbistot(k,ib,iprot),
1088 & " emixsq_pftot",emixsq_pftot(k,ib,iprot),
1089 & " ebis_ftot",ebis_ftot(ib,iprot),
1090 & " emean_ftot",emean_ftot(ib,iprot)
1091 write (iout,*) "k=",k," edecum",edecum,
1092 & " edecum1",edecum1,
1093 & " e2decum",e2decum," ebisdecum",ebisdecum,
1095 & zvdertot(kk,ib,iprot)
1101 ebisdecum=eneps_mix_fbistot(k,ib,iprot)-
1102 & enepsave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1103 edecum=eneps_mix_ftot(k,ib,iprot)
1104 & -enepsave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1105 e2decum=eneps_mixsq_ftot(k,ib,iprot)
1106 & -enepsave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1107 zvepstot(k,ib,iprot)=ww(1)*(Rgas*fac2*(2*edecum
1108 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1109 & *edecum))+fac*ebisdecum)
1111 write (iout,*) ib,k,"ebisdecum",ebisdecum,
1112 & " edecum",edecum," e2decum",e2decum," zvepstot",
1113 & zvepstot(k,ibb,iprot)
1115 & " enepsave_ftot",enepsave_ftot(k,ib,iprot),
1116 & " eneps_mix_ftot",eneps_mix_ftot(k,ib,iprot),
1117 & " eneps_mix_fbistot",eneps_mix_fbistot(k,ib,iprot),
1118 & " eneps_mixsq_ftot",eneps_mixsq_ftot(k,ib,iprot),
1119 & " emean_ftot",emean_ftot(ib,iprot),
1120 & " ebis_ftot",ebis_ftot(ib,iprot)
1124 if (mod_tor .or. mod_sccor) then
1126 ebisdecum=entor_mix_fbistot(k,ib,iprot)-
1127 & entorave_ftot(k,ib,iprot)*ebis_ftot(ib,iprot)
1128 edecum=entor_mix_fprimtot(k,ib,iprot)
1129 & -entorave_fprimtot(k,ib,iprot)*emean_ftot(ib,iprot)
1130 edecum1=entor_mix_ftot(k,ib,iprot)
1131 & -entorave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1132 e2decum=entor_mixsq_ftot(k,ib,iprot)
1133 & -entorave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1134 zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1135 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1136 & *edecum1))-entorave_fbistot(k,ib,iprot)
1140 & " entorave_ftot",entorave_ftot(k,ib,iprot),
1141 & " entor_mix_ftot",entor_mix_ftot(k,ib,iprot),
1142 & " entor_mix_fprimtot",entor_mix_fprimtot(k,ib,iprot),
1143 & " entor_mix_fbistot",entor_mix_fbistot(k,ib,iprot),
1144 & " entorave_fprimtot",entorave_fprimtot(k,ib,iprot),
1145 & " entorave_fbistot",entorave_fbistot(k,ib,iprot),
1146 & " entor_mixsq_ftot",entor_mixsq_ftot(k,ib,iprot),
1147 & " emean_ftot",emean_ftot(ib,iprot),
1148 & " ebis_ftot",ebis_ftot(ib,iprot)
1149 write (iout,*) ib,k," ebisdecum",ebisdecum,
1150 & " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1151 & " zvtortot",zvtortot(k,ibb,iprot)
1159 edecum1=enang_mix_ftot(k,ib,iprot)
1160 & -enangave_ftot(k,ib,iprot)*emean_ftot(ib,iprot)
1161 e2decum=enang_mixsq_ftot(k,ib,iprot)
1162 & -enangave_ftot(k,ib,iprot)*esquare_ftot(ib,iprot)
1163 zvtortot(k,ib,iprot)=Rgas*fac2*(2*edecum
1164 & -fac*(e2decum-2*emean_ftot(ib,iprot)
1165 & *edecum1))-entorave_fbistot(k,ib,iprot)
1169 & " enangave_ftot",enangave_ftot(k,ib,iprot),
1170 & " enang_mix_ftot",enang_mix_ftot(k,ib,iprot),
1171 & " enang_mixsq_ftot",enang_mixsq_ftot(k,ib,iprot),
1172 & " emean_ftot",emean_ftot(ib,iprot)
1173 write (iout,*) ib,k," ebisdecum",ebisdecum,
1174 & " edecum",edecum," edecum1",edecum1," e2decum",e2decum,
1175 & " zvtortot",zvangtot(k,ibb,iprot)
1180 write (iout,*) "iprot",iprot,
1182 write (iout,*) "zvder",
1183 & (zvdertot(k,ibb,iprot),k=1,kk)
1189 c------------------------------------------------------------------------------
1190 subroutine propderiv(ib,inat,iprot)
1191 c Compute the derivatives of conformation-dependent averages in force-field
1194 include "DIMENSIONS"
1195 include "DIMENSIONS.ZSCOPT"
1198 integer IERROR,ERRCODE,status(MPI_STATUS_SIZE)
1201 include "COMMON.WEIGHTS"
1202 include "COMMON.ENERGIES"
1203 include "COMMON.CLASSES"
1204 include "COMMON.IOUNITS"
1205 include "COMMON.DERIV"
1206 include "COMMON.WEIGHTDER"
1207 include "COMMON.VMCPAR"
1208 include "COMMON.GEO"
1209 include "COMMON.OPTIM"
1210 include "COMMON.NAMES"
1211 include "COMMON.COMPAR"
1213 include "COMMON.MPI"
1215 include "COMMON.TIME1"
1216 include "COMMON.VAR"
1217 double precision aux,zscder
1218 double precision fac
1219 double precision efj(2),emj(2),varj(2)
1220 double precision eavpfj(max_ene,2),emixpfj(max_ene,2),
1221 & emixsqpfj(max_ene,2),enepsavefj(nntyp,2),enepsmixfj(nntyp,2),
1222 & enepsmixsqfj(nntyp,2)
1224 integer i,ii,ind,j,k,kk,iprot,ib,ibb,l1
1225 integer iti,itj,jj,icant
1228 integer ind_shield /25/
1229 c write (iout,*) "Entered PROPDERIV"
1231 fac = betaT(ib,inat+2,iprot)
1232 c write (iout,*) "derivatives: iprot",iprot,
1234 do l=1,natdim(inat,iprot)
1237 if (imask(k).gt.0 .and. k.ne.ind_shield) then
1239 rder(kk,l)=-fac*(nu_pf(k,l,ib,inat,iprot)
1240 & -nuave(l,ib,inat,iprot)*
1241 & eave_nat_pftot(k,ib,inat,iprot))
1246 do l=1,natdim(inat,iprot)
1248 reps(k,l)=-ww(1)*fac*(nuepsave_f(k,l,ib,inat,iprot)-
1249 & nuave(l,ib,inat,iprot)*enepsave_nat_ftot(k,ib,inat,iprot))
1253 if (mod_tor .or. mod_sccor) then
1254 do l=1,natdim(inat,iprot)
1256 rtor(k,l)=-fac*(nutorave_f(k,l,ib,inat,iprot)-
1257 & nuave(l,inat,ib,iprot)*entorave_nat_ftot(k,ib,inat,iprot))
1262 do l=1,natdim(inat,iprot)
1264 rang(k,l)=-fac*(nuangave_f(k,l,ib,inat,iprot)-
1265 & nuave(l,inat,ib,iprot)*enangave_nat_ftot(k,ib,inat,iprot))
1271 c-------------------------------------------------------------------------------
1272 subroutine ene_eval(iprot,i,ii)
1274 include "DIMENSIONS"
1275 include "DIMENSIONS.ZSCOPT"
1278 integer IERROR,ErrCode
1279 include "COMMON.MPI"
1281 include "COMMON.CONTROL"
1282 include "COMMON.COMPAR"
1283 include "COMMON.WEIGHTS"
1284 include "COMMON.WEIGHTDER"
1285 include "COMMON.CLASSES"
1286 include "COMMON.ENERGIES"
1287 include "COMMON.IOUNITS"
1288 include "COMMON.VMCPAR"
1289 include "COMMON.NAMES"
1290 include "COMMON.INTERACT"
1291 include "COMMON.TIME1"
1292 include "COMMON.CHAIN"
1293 include "COMMON.VAR"
1294 include "COMMON.GEO"
1295 include "COMMON.SCCOR"
1296 include "COMMON.TORSION"
1297 include "COMMON.LOCAL"
1298 include "COMMON.FFIELD"
1299 C Define local variables
1300 integer i,ii,ij,jj,j,k,ik,jk,l,ll,lll,iprot,iblock
1301 double precision energia(0:max_ene)
1302 double precision etoti,enepsjk
1303 double precision ftune_eps
1308 integer ind_shield /25/
1312 write (iout,*) "ene_eval: iprot",iprot," i",i," ii",ii
1314 if (init_ene .or. (mod_fourier(nloctyp).gt.0
1315 & .or. mod_elec .or. mod_scp .or. mod_tor .or. mod_sccor
1316 & .or. mod_ang .or. imask(ind_shield).gt.0) ) then
1318 write (iout,*) "FUNC1: doing UNRES"
1319 write (iout,*) "Protein",iprot," i",i," ii",ii," nres",nres
1322 call restore_ccoords(iprot,ii)
1323 call int_from_cart1(.false.)
1325 write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1326 write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1327 write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1328 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1329 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1330 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1331 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1332 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1337 call etotal(energia(0))
1338 c if (natlike(iprot).gt.0)
1339 c & call natsimil(i,ii,iprot,.false.)
1340 e_total(i,iprot)=energia(0)
1342 enetb(ii,j,iprot)=energia(j)
1346 write (iout,*) "ww",ww(:n_ene)
1347 write (iout,'(2i5,18(1pe12.4))') i,ii,
1348 & (enetb(ii,j,iprot),j=1,n_ene)
1349 write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1350 & i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)
1351 call enerprint(energia(0))
1355 write (iout,*) "FUN nntyp",nntyp," eneps",eneps_temp(1,45)
1358 eneps(k,j,ii,iprot)=eneps_temp(k,j)
1362 c write (iout,*) "Matrix eneps"
1365 c write (iout,*) j,k,icant(j,k),
1366 c & eneps_temp(1,icant(j,k)),
1367 c & eneps_temp(2,icant(j,k))
1373 if (mod_tor .or. mod_sccor) then
1374 c write (iout,*) "etor",enetb(ii,13,iprot),enetb(ii,19,iprot)
1377 do ij=-ntortyp+1,ntortyp-1
1378 do jj=-ntortyp+1,ntortyp-1
1379 if (mask_tor(0,jj,ij,iblock).eq.1) then
1381 entor(k,ii,iprot)=etor_temp(0,0,jj,ij,iblock)
1382 c write (iout,*) "etor_temp",restyp(itortyp(j)),
1383 c & restyp(itortyp(i)),
1384 c & k,etor_temp(ji,ii)
1385 else if (mask_tor(0,jj,ij,iblock).gt.1) then
1386 if (tor_mode.eq.0) then
1387 do l=1,nterm(jj,ij,iblock)
1389 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1391 if (mask_tor(0,jj,ij,iblock).gt.2) then
1392 do l=nterm(jj,ij,iblock)+1,2*nterm(jj,ij,iblock)
1394 entor(k,ii,iprot)=etor_temp(l,0,jj,ij,iblock)
1397 else if (tor_mode.eq.1) then
1399 do l=1,nterm_kcc(jj,ij)
1400 do ll=1,nterm_kcc_Tb(jj,ij)
1401 do lll=1,nterm_kcc_Tb(jj,ij)
1403 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1407 if (mask_tor(0,jj,ij,1).gt.2) then
1408 do l=nterm_kcc(jj,ij)+1,2*nterm_kcc(jj,ij)
1409 do ll=1,nterm_kcc_Tb(jj,ij)
1410 do lll=1,nterm_kcc_Tb(jj,ij)
1412 entor(k,ii,iprot)=etor_temp_kcc(lll,ll,l,ij,jj)
1423 do ij=-nsccortyp+1,nsccortyp-1
1424 do jj=-nsccortyp+1,nsccortyp-1
1426 if (mask_tor(ll,jj,ij,1).eq.1) then
1428 entor(k,ii,iprot)=etor_temp(0,ll,jj,ij,1)
1429 c write (iout,*) "etor_temp",restyp(isccortyp(jj)),
1430 c & restyp(isccortyp(ij)),
1431 c & k,etor_temp(jj,ij)
1432 else if (mask_tor(ll,jj,ij,1).gt.1) then
1433 do l=1,nterm_sccor(jj,ij)
1435 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1437 do l=nterm_sccor(jj,ij)+1,2*nterm_sccor(jj,ij)
1439 entor(k,ii,iprot)=etor_temp(l,ll,jj,ij,1)
1446 c AL 8/7/17: optimization of the bending potentials
1450 if (mask_ang(ij).eq.0) cycle
1451 do jj=1,nbend_kcc_Tb(ij)
1453 enbend(k,ii,iprot)=ebend_temp_kcc(jj,ij)*wang
1454 c write (iout,*) "ij",ij," jj",jj," k",k," ebend_temp_kcc",
1455 c & ebend_temp_kcc(jj,ij)
1461 write (iout,*) "Conformation",i
1462 c write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1463 c write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1464 c write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1465 c write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1466 c write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1467 c write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1468 c write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1469 c write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1470 call enerprint(energia(0))
1471 c write(iout,'(4hENER,i10,3f15.3,2i10,i5)')
1472 c & i,energia(0),eini(i,iprot),entfac(i,iprot)
1473 c write (iout,'(i5,18(1pe12.4))') ii,
1474 c & (enetb(ii,j,iprot),j=1,n_ene)
1475 c write (iout,'(i5,16(1pe12.4))') i,
1476 c & (energia(j),j=1,n_ene),energia(0)
1477 c print *,"Processor",me,me1," computing energies",i
1480 if (energia(0).ge.1.0d99) then
1481 write (iout,*) "FUNC1:CHUJ NASTAPIL in energy evaluation for",
1482 & " point",i,". Probably NaNs in some of the energy components."
1483 write (iout,*) "The components of the energy are:"
1484 call enerprint(energia(0))
1485 write (iout,*) "Conformation",i,ii
1486 write (iout,'(8f10.5)') ((c(j,k),j=1,3),k=1,nres)
1487 write (iout,'(8f10.5)') ((c(j,k+nres),j=1,3),k=nnt,nct)
1488 write (iout,'(8f10.4)') (vbld(k),k=2,nres)
1489 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1490 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1491 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1492 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1493 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1494 write (iout,*) "Calculation terminated at this point.",
1495 & " Check the database of conformations"
1498 call MPI_Abort(MPI_COMM_WORLD,ErrCode,IERROR)
1500 stop "SEVERE error in energy calculation"
1504 enetb(ii,1,iprot)=0.0d0
1507 enetb(ii,1,iprot)=enetb(ii,1,iprot)+
1508 & ftune_eps(eps(j,k))*eneps(1,icant(j,k),ii,iprot)+
1509 & eps(j,k)*eneps(2,icant(j,k),ii,iprot)
1515 etoti=etoti+enetb(ii,j,iprot)*ww(j)
1519 write (iout,*) "i",i," ww",ww
1520 write (iout,*) "enetb",(enetb(ii,j,iprot),j=1,n_ene)
1523 e_total(i,iprot)=etoti
1524 c write (iout,*) "FUNC1 natlike",iprot,natlike(iprot)
1525 if (natlike(iprot).gt.0) then
1526 call restore_ccoords(iprot,ii)
1527 call int_from_cart1(.false.)
1528 c call natsimil(i,ii,iprot,.false.)
1530 c & "natsimil",i,ii,(nu(k,ii,iprot),k=1,natlike(iprot))
1533 write (iout,'(i5,18(1pe12.4))') i,
1534 & (enetb(ii,j,iprot),j=1,n_ene)
1535 write (iout,'(18(1pe12.4))')
1536 & (eneps(1,j,ii,iprot),j=201,210)
1537 write(iout,'(4hENER,i10,3f15.3)')
1538 & i,e_total(i,iprot),eini(i,iprot),entfac(i,iprot)