source/wham/src-HCD-5D/enecalc1.F

   1       subroutine enecalc(islice,*)
   2       implicit none
   3       include "DIMENSIONS"
   4       include "DIMENSIONS.ZSCOPT"
   5       include "DIMENSIONS.FREE"
   6 #ifdef MPI
   7       include "mpif.h"
   8       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
   9       include "COMMON.MPI"
  10 #endif
  11       include "COMMON.CHAIN"
  12       include "COMMON.IOUNITS"
  13       include "COMMON.PROTFILES"
  14       include "COMMON.NAMES"
  15       include "COMMON.VAR"
  16       include "COMMON.SBRIDGE"
  17       include "COMMON.GEO"
  18       include "COMMON.FFIELD"
  19       include "COMMON.ENEPS"
  20       include "COMMON.LOCAL"
  21       include "COMMON.WEIGHTS"
  22       include "COMMON.INTERACT"
  23       include "COMMON.FREE"
  24       include "COMMON.ENERGIES"
  25       include "COMMON.CONTROL"
  26       include "COMMON.TORCNSTR"
  27       character*64 nazwa
  28       character*80 bxname
  29       character*3 liczba
  30       double precision qwolynes
  31       external qwolynes
  32       integer errmsg_count,maxerrmsg_count /100/
  33       double precision rmsnat,gyrate
  34       external rmsnat,gyrate
  35 c      double precision tole /1.0d-1/
  36       integer i,itj,ii,iii,j,k,l,licz,ipermin
  37       integer ir,ib,ipar,iparm
  38       integer iscor,islice,scount_buff(0:99)
  39       real*4 csingle(3,maxres2)
  40       double precision energ
  41       double precision temp
  42       integer ilen,iroof
  43       external ilen,iroof
  44       double precision energia(0:max_ene),rmsdev,efree,eini
  45       double precision fT(6),quot,quotl,kfacl,kfac /2.4d0/,T0 /3.0d2/
  46       double precision tt
  47       integer snk_p(MaxR,MaxT_h,Max_parm)
  48       logical lerr
  49       integer ncont,icont(2,maxcont),isecstr(maxres)
  50       character*256 bprotfile_temp
  51       double precision totlength
  52       call opentmp(islice,ientout,bprotfile_temp)
  53       iii=0
  54       ii=0
  55       errmsg_count=0
  56 c      write (iout,*) "enecalc: nparmset ",nparmset
  57 c      write (iout,*) "enecalc: tormode ",tor_mode
  58       write (iout,*) "ns",ns," dyn_ss",dyn_ss,(iss(i),i=1,ns)
  59       if (ns.gt.0.and.dyn_ss) then
  60           do i=nss+1,nhpb
  61             ihpb(i-nss)=ihpb(i)
  62             jhpb(i-nss)=jhpb(i)
  63             forcon(i-nss)=forcon(i)
  64             dhpb(i-nss)=dhpb(i)
  65           enddo
  66           nhpb=nhpb-nss
  67           nss=0
  68           call hpb_partition
  69           do i=1,ns
  70             dyn_ss_mask(iss(i))=.true.
  71           enddo
  72       endif
  73       write (iout,*) "dyn_ss_mask",(dyn_ss_mask(i),i=1,nres)
  74 #ifdef MPI
  75       do iparm=1,nParmSet
  76         do ib=1,nT_h(iparm)
  77           do i=1,nR(ib,iparm)
  78             snk_p(i,ib,iparm)=0
  79           enddo
  80         enddo
  81       enddo
  82       write (iout,*) "indstart(me1),indend(me1)"
  83      &,indstart(me1),indend(me1)
  84       do i=indstart(me1),indend(me1)
  85 #else
  86       do iparm=1,nParmSet
  87         do ib=1,nT_h(iparm)
  88           do i=1,nR(ib,iparm)
  89             snk(i,ib,iparm)=0
  90           enddo
  91         enddo
  92       enddo
  93       do i=1,ntot
  94 #endif
  95
  96         read(ientout,rec=i,err=101)
  97      &    ((csingle(l,k),l=1,3),k=1,nres),
  98      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
  99      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 100      &    eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
 101          if (indpdb.gt.0) then
 102            do k=1,nres
 103              do l=1,3
 104                c(l,k)=csingle(l,k)
 105              enddo
 106            enddo
 107            do k=nnt,nct
 108              do l=1,3
 109                c(l,k+nres)=csingle(l,k+nres)
 110              enddo
 111            enddo
 112            anatemp= 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 113            q(nQ+1,iii+1)=rmsnat(iii+1,ipermin)
 114          endif
 115 c         write (iout,*) iii+1,q(nQ+3,iii+1),q(nQ+4,iii+1),q(nQ+5,iii+1)
 116 c        fT=T0*beta_h(ib,ipar)*1.987D-3
 117 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
 118         if (rescale_mode.eq.1) then
 119           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 120 #if defined(FUNCTH)
 121           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 122           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 123 #elif defined(FUNCT)
 124           ft(6)=quot
 125 #else
 126           ft(6)=1.0d0
 127 #endif
 128           quotl=1.0d0
 129           kfacl=1.0d0
 130           do l=1,5
 131             quotl=quotl*quot
 132             kfacl=kfacl*kfac
 133             fT(l)=kfacl/(kfacl-1.0d0+quotl)
 134           enddo
 135         else if (rescale_mode.eq.2) then
 136           quot=1.0d0/(T0*beta_h(ib,ipar)*1.987D-3)
 137 #if defined(FUNCTH)
 138           tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
 139           ft(6)=(320.0+80.0*dtanh((tt-320.0)/80.0))/320.0
 140 #elif defined(FUNCT)
 141           ft(6)=quot
 142 #else
 143           ft(6)=1.0d0
 144 #endif
 145           quotl=1.0d0
 146           do l=1,5
 147             quotl=quotl*quot
 148             fT(l)=1.12692801104297249644d0/
 149      &         dlog(dexp(quotl)+dexp(-quotl))
 150           enddo
 151         else if (rescale_mode.eq.0) then
 152           do l=1,5
 153             fT(l)=1.0d0
 154           enddo
 155         else
 156           write (iout,*) "Error in ECECALC: wrong RESCALE_MODE",
 157      &     rescale_mode
 158           call flush(iout)
 159           return1
 160         endif
 161
 162 c        write (iout,*) "T",1.0d0/(beta_h(ib,ipar)*1.987D-3)," T0",T0,
 163 c     &   " kfac",kfac,"quot",quot," fT",fT
 164         do j=1,2*nres
 165           do k=1,3
 166             c(k,j)=csingle(k,j)
 167           enddo
 168         enddo
 169         call int_from_cart1(.false.)
 170         ii=ii+1
 171         do iparm=1,nparmset
 172
 173         call restore_parm(iparm)
 174 #ifdef DEBUG
 175             write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 176      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 177      &      wtor_d,wsccor,wbond
 178 #endif
 179 c        write (iout,*) "tuz przed energia"
 180         call etotal(energia(0),fT)
 181 c        write (iout,*) "tuz za energia"
 182 #ifdef DEBUG
 183         write (iout,*) "Conformation",i
 184         write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
 185      &                            ((c(l,k+nres),l=1,3),k=nnt,nct)
 186         call enerprint(energia(0),fT)
 187 c        write (iout,'(2i5,21f8.2)') i,iparm,(energia(k),k=1,21)
 188 c        write (iout,*) "ftors(1)",ftors(1)
 189 c        call briefout(i,energia(0))
 190 c        temp=1.0d0/(beta_h(ib,ipar)*1.987D-3)
 191 c        write (iout,*) "temp", temp
 192 c        call pdbout(i,temp,energia(0),energia(0),0.0d0,0.0d0)
 193 #endif
 194         if (isnan(energia(0)) .or. energia(1).ge.1.0d20
 195      &     .or. energia(0).ge.1.0d20) then
 196           write (iout,*) "NaNs detected in some of the energy",
 197      &     " components for conformation",ii+1
 198           write (iout,*) "The Cartesian geometry is:"
 199           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 200           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 201           write (iout,*) "The internal geometry is:"
 202 c          call intout
 203 c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 204           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 205           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 206           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 207           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 208           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 209           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 210           write (iout,*) "The components of the energy are:"
 211           call enerprint(energia(0),fT)
 212           write (iout,*)
 213      &      "This conformation WILL NOT be added to the database."
 214           call flush(iout)
 215           goto 121
 216         else
 217 #ifdef DEBUG
 218           if (ipar.eq.iparm) write (iout,*) i,iparm,
 219      &      1.0d0/(beta_h(ib,ipar)*1.987D-3),eini,energia(0)
 220 #endif
 221 c          write (iout,*) "eini",eini,"energia(0)",energia(0)," diff",
 222 c     &       eini-energia(0)
 223           if (ipar.eq.iparm .and. einicheck.gt.0 .and.
 224 !     &      dabs(eini-energia(0)-energia(27)).gt.tole) then
 225      &      dabs(eini-energia(0)).gt.tole) then
 226             if (errmsg_count.le.maxerrmsg_count) then
 227               write (iout,'(2a,2e15.5,a,2i8,a,f8.1)')
 228      &         "Warning: energy differs remarkably from ",
 229 !     &      " the value read in: ",energia(0)+energia(27),eini," point",
 230      &      " the value read in: ",energia(0),eini," point",
 231      &         iii+1,indstart(me1)+iii," T",
 232      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
 233           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
 234      &                            ((c(l,k+nres),l=1,3),k=nnt,nct)
 235 c              call intout
 236 c              call pdbout(indstart(me1)+iii,
 237 c     & 1.0d0/(1.987D-3*beta_h(ib,ipar)),energia(0),eini,0.0d0,0.0d0)
 238               call enerprint(energia(0),fT)
 239               errmsg_count=errmsg_count+1
 240               if (errmsg_count.gt.maxerrmsg_count)
 241      &          write (iout,*) "Too many warning messages"
 242               if (einicheck.gt.1) then
 243                 write (iout,*) "Calculation stopped."
 244                 call flush(iout)
 245 #ifdef MPI
 246                 call MPI_Abort(WHAM_COMM,IERROR,ERRCODE)
 247 #endif
 248                 call flush(iout)
 249                 return1
 250               endif
 251             endif
 252           endif
 253 C          write (iout,*) "Czy tu dochodze"
 254           potE(iii+1,iparm)=energia(0)
 255           do k=1,n_ene
 256             enetb(k,iii+1,iparm)=energia(k)
 257           enddo
 258 #ifdef DEBUG
 259           write (iout,'(2i5,f10.1,3e15.5)') i,iii,
 260      &     1.0d0/(beta_h(ib,ipar)*1.987D-3),energia(0),eini,efree
 261 c          call enerprint(energia(0),fT)
 262 #endif
 263 #ifdef DEBUG
 264           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 265           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 266           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 267           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 268           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 269           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 270           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 271           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 272           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 273           write (iout,'(8f10.5)') (q(k,iii+1),k=1,nQ)
 274           write (iout,'(f10.5,i10)') rmsdev,iscor
 275           call enerprint(energia(0),fT)
 276         call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
 277 #endif
 278         endif
 279
 280         enddo ! iparm
 281         q(nQ+2,iii+1)=gyrate(iii+1)
 282 c 8/28/2020 Adam - determine the fraction of secondary structures.
 283         call elecont(.false.,ncont,icont,nnt,nct-1,1)
 284         call secondary2(.false.,.false.,ncont,icont,isecstr)
 285 #ifdef DEBUG
 286         write (iout,*) "secondary structure"
 287         write (iout,'(80i1)') (isecstr(k),k=1,nres)
 288 #endif
 289         q(nQ+3,iii+1)=0.0d0
 290         q(nQ+4,iii+1)=0.0d0
 291         q(nQ+5,iii+1)=0.0d0
 292         totlength=0.0d0
 293         do k=nnt,nct
 294           if (itype(k).eq.ntyp1) cycle
 295           totlength=totlength+1.0d0
 296           l=isecstr(k)
 297           q(nQ+3+l,iii+1)=q(nQ+3+l,iii+1)+1.0d0
 298         enddo
 299         q(nQ+3,iii+1)=q(nQ+3,iii+1)/totlength
 300         q(nQ+4,iii+1)=q(nQ+4,iii+1)/totlength
 301         q(nQ+5,iii+1)=q(nQ+5,iii+1)/totlength
 302 c        write (iout,*) "iii",iii," nssbond",nssbond,nss
 303 c        q(nQ+6,iii+1)=nssbond
 304         q(nQ+6,iii+1)=nss
 305
 306         iii=iii+1
 307         if (q(1,iii).le.0.0d0 .and. indpdb.gt.0)
 308      &    q(1,iii)=qwolynes(0,0,ipermin)
 309 c        write (iout,*) "iii",iii," q",q(1,iii)
 310         write (ientout,rec=iii)
 311      &   ((csingle(l,k),l=1,3),k=1,nres),
 312      &   ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 313      &   nss,(ihpb(k),jhpb(k),k=1,nss),
 314      &   potE(iii,ipar),efree,rmsdev,(q(k,iii),k=1,nQ),iR,ib,ipar
 315 c        write (iout,'(2i5,2e15.5)') ii,iii,potE(iii,ipar),efree
 316 #ifdef MPI
 317         if (separate_parset) then
 318           snk_p(iR,ib,1)=snk_p(iR,ib,1)+1
 319         else
 320           snk_p(iR,ib,ipar)=snk_p(iR,ib,ipar)+1
 321         endif
 322 c        write (iout,*) "iii",iii," iR",iR," ib",ib," ipar",ipar,
 323 c     &   " snk",snk_p(iR,ib,ipar)
 324 #else
 325         snk(iR,ib,ipar,islice)=snk(iR,ib,ipar,islice)+1
 326 #endif
 327   121   continue
 328       enddo
 329 #ifdef MPI
 330       scount_buff(me)=iii
 331       write (iout,*) "Me",me," scount_buff",scount_buff(me)
 332       call flush(iout)
 333 c  Master gathers updated numbers of conformations written by all procs.
 334 c      call MPI_AllGather(MPI_IN_PLACE,1,MPI_DATATYPE_NULL,scount(0),1,
 335 c     &  MPI_INTEGER, WHAM_COMM, IERROR)
 336       call MPI_AllGather( scount_buff(me), 1, MPI_INTEGER, scount(0), 1,
 337      &  MPI_INTEGER, WHAM_COMM, IERROR)
 338
 339       indstart(0)=1
 340       indend(0)=scount(0)
 341       do i=1, Nprocs-1
 342         indstart(i)=indend(i-1)+1
 343         indend(i)=indstart(i)+scount(i)-1
 344       enddo
 345       write (iout,*)
 346       write (iout,*) "Revised conformation counts"
 347       do i=0,nprocs1-1
 348         write (iout,'(a,i5,a,i7,a,i7,a,i7)')
 349      &    "Processor",i," indstart",indstart(i),
 350      &    " indend",indend(i)," count",scount(i)
 351       enddo
 352       call flush(iout)
 353       call MPI_AllReduce(snk_p(1,1,1),snk(1,1,1,islice),
 354      &  MaxR*MaxT_h*nParmSet,
 355      &  MPI_INTEGER,MPI_SUM,WHAM_COMM,IERROR)
 356 #endif
 357       stot(islice)=0
 358       do iparm=1,nParmSet
 359         do ib=1,nT_h(iparm)
 360           do i=1,nR(ib,iparm)
 361             stot(islice)=stot(islice)+snk(i,ib,iparm,islice)
 362           enddo
 363         enddo
 364       enddo
 365       write (iout,*) "Revised SNK"
 366       do iparm=1,nParmSet
 367         do ib=1,nT_h(iparm)
 368           write (iout,'("Param",i3," Temp",f6.1,3x,32i8)')
 369      &     iparm,1.0d0/(1.987D-3*beta_h(ib,iparm)),
 370      &     (snk(i,ib,iparm,islice),i=1,nR(ib,iparm))
 371           write (iout,*) "snk_p",(snk_p(i,ib,iparm),i=1,nR(ib,iparm))
 372         enddo
 373       enddo
 374       write (iout,'("Total",i10)') stot(islice)
 375       call flush(iout)
 376       return
 377   101 write (iout,*) "Error in scratchfile."
 378       call flush(iout)
 379       return1
 380       end
 381 c------------------------------------------------------------------------------
 382       subroutine write_dbase(islice,*)
 383       implicit none
 384       include "DIMENSIONS"
 385       include "DIMENSIONS.ZSCOPT"
 386       include "DIMENSIONS.FREE"
 387       include "DIMENSIONS.COMPAR"
 388 #ifdef MPI
 389       include "mpif.h"
 390       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 391       include "COMMON.MPI"
 392 #endif
 393       include "COMMON.CONTROL"
 394       include "COMMON.CHAIN"
 395       include "COMMON.IOUNITS"
 396       include "COMMON.PROTFILES"
 397       include "COMMON.NAMES"
 398       include "COMMON.VAR"
 399       include "COMMON.SBRIDGE"
 400       include "COMMON.GEO"
 401       include "COMMON.FFIELD"
 402       include "COMMON.ENEPS"
 403       include "COMMON.LOCAL"
 404       include "COMMON.WEIGHTS"
 405       include "COMMON.INTERACT"
 406       include "COMMON.FREE"
 407       include "COMMON.ENERGIES"
 408       include "COMMON.COMPAR"
 409       include "COMMON.PROT"
 410       include "COMMON.CONTACTS1"
 411       character*64 nazwa
 412       character*80 bxname,cxname
 413       character*256 bprotfile_temp
 414       character*3 liczba,licz
 415       character*2 licz2
 416       integer i,itj,ii,iii,j,k,l
 417       integer ixdrf,iret
 418       integer iscor,islice
 419       double precision rmsdev,efree,eini
 420       real*4 csingle(3,maxres2)
 421       double precision energ
 422       integer ilen,iroof
 423       external ilen,iroof
 424       integer ir,ib,iparm, scount_buff(0:99)
 425       integer isecstr(maxres)
 426       write (licz2,'(bz,i2.2)') islice
 427       call opentmp(islice,ientout,bprotfile_temp)
 428       write (iout,*) "bprotfile_temp ",bprotfile_temp
 429       call flush(iout)
 430       if (.not.bxfile .and. .not. cxfile .and. indpdb.eq.0
 431      &   .and. ensembles.eq.0) then
 432         close(ientout,status="delete")
 433         return
 434       endif
 435 #ifdef MPI
 436       write (liczba,'(bz,i3.3)') me
 437       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 438         if (.not.separate_parset) then
 439           bxname = prefix(:ilen(prefix))//liczba//".bx"
 440         else
 441           write (licz,'(bz,i3.3)') myparm
 442           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 443         endif
 444         open (ientin,file=bxname,status="unknown",
 445      &    form="unformatted",access="direct",recl=lenrec1)
 446       endif
 447 #else
 448       if (bxfile .or. cxfile .or. ensembles.gt.0) then
 449         if (nslice.eq.1) then
 450           bxname = prefix(:ilen(prefix))//".bx"
 451         else
 452           bxname = prefix(:ilen(prefix))//
 453      &           "_slice_"//licz2//".bx"
 454         endif
 455         open (ientin,file=bxname,status="unknown",
 456      &    form="unformatted",access="direct",recl=lenrec1)
 457         write (iout,*) "Calculating energies; writing geometry",
 458      & " and energy components to ",bxname(:ilen(bxname))
 459       endif
 460 #if (defined(AIX) && !defined(JUBL))
 461         call xdrfopen_(ixdrf,cxname, "w", iret)
 462 #else
 463         call xdrfopen(ixdrf,cxname, "w", iret)
 464 #endif
 465         if (iret.eq.0) then
 466           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 467           cxfile=.fale.
 468         endif
 469       endif
 470 #endif
 471       if (indpdb.gt.0) then
 472         if (nslice.eq.1) then
 473 #ifdef MPI
 474          if (.not.separate_parset) then
 475            statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 476      &       //liczba//'.stat'
 477          else
 478            write (licz,'(bz,i3.3)') myparm
 479            statname=prefix(:ilen(prefix))//'_par'//licz//'_'//
 480      &      pot(:ilen(pot))//liczba//'.stat'
 481          endif
 482
 483 #else
 484           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 485 #endif
 486         else
 487 #ifdef MPI
 488          if (.not.separate_parset) then
 489           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 490      &      "_slice_"//licz2//liczba//'.stat'
 491          else
 492           write (licz,'(bz,i3.3)') myparm
 493           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//
 494      &      '_par'//licz//"_slice_"//licz2//liczba//'.stat'
 495          endif
 496 #else
 497           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))
 498      &      //"_slice_"//licz2//'.stat'
 499 #endif
 500         endif
 501         open(istat,file=statname,status="unknown")
 502       endif
 503
 504 #ifdef MPI
 505       do i=1,scount(me)
 506 #else
 507       do i=1,ntot(islice)
 508 #endif
 509         read(ientout,rec=i,err=101)
 510      &    ((csingle(l,k),l=1,3),k=1,nres),
 511      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 512      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 513      &    eini,efree,rmsdev,(q(k,i),k=1,nQ),iR,ib,iparm
 514 c        write (iout,*) iR,ib,iparm,eini,efree
 515         do j=1,2*nres
 516           do k=1,3
 517             c(k,j)=csingle(k,j)
 518           enddo
 519         enddo
 520         call int_from_cart1(.false.)
 521         iscore=0
 522 c        write (iout,*) "Calling conf_compar",i
 523 c        call flush(iout)
 524          anatemp= 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 525         if (indpdb.gt.0) then
 526           call conf_compar(i,.false.,.true.)
 527 c        else
 528 c            call elecont(.false.,ncont,icont,nnt,nct)
 529 c            call secondary2(.false.,.false.,ncont,icont,isecstr)
 530         endif
 531 c        write (iout,*) "Exit conf_compar",i
 532 c        call flush(iout)
 533         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i)
 534      &    ((csingle(l,k),l=1,3),k=1,nres),
 535      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 536      &    nss,(ihpb(k),jhpb(k),k=1,nss),
 537 c     &    potE(i,iparm),-entfac(i),rms_nat,iscore
 538      &    potE(i,nparmset),-entfac(i),rms_nat,iscore
 539 c        write (iout,'(2i5,3e15.5)') i,me,potE(i,iparm),-entfac(i)
 540 #ifndef MPI
 541         if (cxfile) call cxwrite(ixdrf,csingle,potE(i,nparmset),
 542      &    -entfac(i),rms_nat,iscore)
 543 #endif
 544       enddo
 545       close(ientout,status="delete")
 546       close(istat)
 547       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 548 #ifdef MPI
 549       call MPI_Barrier(WHAM_COMM,IERROR)
 550       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile
 551      &   .and. ensembles.eq.0) return
 552       write (iout,*)
 553       if (bxfile .or. ensembles.gt.0) then
 554         if (nslice.eq.1) then
 555           if (.not.separate_parset) then
 556             bxname = prefix(:ilen(prefix))//".bx"
 557           else
 558             write (licz,'(bz,i3.3)') myparm
 559             bxname = prefix(:ilen(prefix))//"_par"//licz//".bx"
 560           endif
 561         else
 562           if (.not.separate_parset) then
 563             bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
 564           else
 565             write (licz,'(bz,i3.3)') myparm
 566             bxname = prefix(:ilen(prefix))//"par_"//licz//
 567      &        "_slice_"//licz2//".bx"
 568           endif
 569         endif
 570         open (ientout,file=bxname,status="unknown",
 571      &      form="unformatted",access="direct",recl=lenrec1)
 572         write (iout,*) "Master is creating binary database ",
 573      &   bxname(:ilen(bxname))
 574       endif
 575       if (cxfile) then
 576         if (nslice.eq.1) then
 577           if (.not.separate_parset) then
 578             cxname = prefix(:ilen(prefix))//".cx"
 579           else
 580             cxname = prefix(:ilen(prefix))//"_par"//licz//".cx"
 581           endif
 582         else
 583           if (.not.separate_parset) then
 584             cxname = prefix(:ilen(prefix))//
 585      &             "_slice_"//licz2//".cx"
 586           else
 587             cxname = prefix(:ilen(prefix))//"_par"//licz//
 588      &             "_slice_"//licz2//".cx"
 589           endif
 590         endif
 591 #if (defined(AIX) && !defined(JUBL))
 592         call xdrfopen_(ixdrf,cxname, "w", iret)
 593 #else
 594         call xdrfopen(ixdrf,cxname, "w", iret)
 595 #endif
 596         if (iret.eq.0) then
 597           write (iout,*) "Error opening cxfile ",cxname(:ilen(cxname))
 598           cxfile=.false.
 599         endif
 600       endif
 601       do j=0,nprocs-1
 602         write (liczba,'(bz,i3.3)') j
 603         if (separate_parset) then
 604           write (licz,'(bz,i3.3)') myparm
 605           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
 606         else
 607           bxname = prefix(:ilen(prefix))//liczba//".bx"
 608         endif
 609         open (ientin,file=bxname,status="unknown",
 610      &    form="unformatted",access="direct",recl=lenrec1)
 611         write (iout,*) "Master is reading conformations from ",
 612      &   bxname(:ilen(bxname))
 613         iii = 0
 614 c        write (iout,*) "j",j," indstart",indstart(j)," indend",indend(j)
 615 c        call flush(iout)
 616         do i=indstart(j),indend(j)
 617           iii = iii+1
 618           read(ientin,rec=iii,err=101)
 619      &      ((csingle(l,k),l=1,3),k=1,nres),
 620      &      ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 621      &      nss,(ihpb(k),jhpb(k),k=1,nss),
 622      &      eini,efree,rmsdev,iscor
 623           if (bxfile .or. ensembles.gt.0) then
 624             write (ientout,rec=i)
 625      &        ((csingle(l,k),l=1,3),k=1,nres),
 626      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 627      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 628      &        eini,efree,rmsdev,iscor
 629           endif
 630           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 631 #ifdef DEBUG
 632           do k=1,2*nres
 633             do l=1,3
 634               c(l,k)=csingle(l,k)
 635             enddo
 636           enddo
 637           call int_from_cart1(.false.)
 638           write (iout,'(2i5,3e15.5)') i,iii,eini,efree
 639           write (iout,*) "The Cartesian geometry is:"
 640           write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 641           write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 642           write (iout,*) "The internal geometry is:"
 643           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 644           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 645           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 646           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 647           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 648           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 649           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 650           write (iout,'(f10.5,i5)') rmsdev,iscor
 651 #endif
 652         enddo ! i
 653         write (iout,*) iii," conformations (from",indstart(j)," to",
 654      &   indend(j),") read from ",
 655      &   bxname(:ilen(bxname))
 656         close (ientin,status="delete")
 657       enddo ! j
 658       if (bxfile .or. cxfile .or. ensembles.gt.0) close (ientout)
 659 #if (defined(AIX) && !defined(JUBL))
 660       if (cxfile) call xdrfclose_(ixdrf,cxname,iret)
 661 #else
 662       if (cxfile) call xdrfclose(ixdrf,cxname,iret)
 663 #endif
 664 #endif
 665       return
 666   101 write (iout,*) "Error in scratchfile."
 667       call flush(iout)
 668       return1
 669       end
 670 c-------------------------------------------------------------------------------
 671       subroutine cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
 672       implicit none
 673       include "DIMENSIONS"
 674       include "DIMENSIONS.ZSCOPT"
 675       include "DIMENSIONS.FREE"
 676       include "DIMENSIONS.COMPAR"
 677 #ifdef MPI
 678       include "mpif.h"
 679       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 680       include "COMMON.MPI"
 681 #endif
 682       include "COMMON.CONTROL"
 683       include "COMMON.CHAIN"
 684       include "COMMON.IOUNITS"
 685       include "COMMON.PROTFILES"
 686       include "COMMON.NAMES"
 687       include "COMMON.VAR"
 688       include "COMMON.SBRIDGE"
 689       include "COMMON.GEO"
 690       include "COMMON.FFIELD"
 691       include "COMMON.ENEPS"
 692       include "COMMON.LOCAL"
 693       include "COMMON.WEIGHTS"
 694       include "COMMON.INTERACT"
 695       include "COMMON.FREE"
 696       include "COMMON.ENERGIES"
 697       include "COMMON.COMPAR"
 698       include "COMMON.PROT"
 699       integer i,j,itmp,iscor,iret,ixdrf
 700       double precision rmsdev,efree,eini
 701       real*4 csingle(3,maxres2),xoord(3,maxres2+2)
 702       real*4 prec
 703
 704 c      write (iout,*) "cxwrite"
 705 c      call flush(iout)
 706       prec=10000.0
 707       do i=1,nres
 708        do j=1,3
 709         xoord(j,i)=csingle(j,i)
 710        enddo
 711       enddo
 712       do i=nnt,nct
 713        do j=1,3
 714         xoord(j,nres+i-nnt+1)=csingle(j,i+nres)
 715        enddo
 716       enddo
 717
 718       itmp=nres+nct-nnt+1
 719
 720 c      write (iout,*) "itmp",itmp
 721 c      call flush(iout)
 722 #if (defined(AIX) && !defined(JUBL))
 723       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 724
 725 c      write (iout,*) "xdrf3dfcoord"
 726 c      call flush(iout)
 727       call xdrfint_(ixdrf, nss, iret)
 728       do j=1,nss
 729            if (dyn_ss) then
 730             call xdrfint(ixdrf, idssb(j), iret)
 731             call xdrfint(ixdrf, jdssb(j), iret)
 732            else
 733             call xdrfint_(ixdrf, ihpb(j), iret)
 734             call xdrfint_(ixdrf, jhpb(j), iret)
 735            endif
 736       enddo
 737       call xdrffloat_(ixdrf,real(eini),iret)
 738       call xdrffloat_(ixdrf,real(efree),iret)
 739       call xdrffloat_(ixdrf,real(rmsdev),iret)
 740       call xdrfint_(ixdrf,iscor,iret)
 741 #else
 742       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 743
 744       call xdrfint(ixdrf, nss, iret)
 745       do j=1,nss
 746            if (dyn_ss) then
 747             call xdrfint(ixdrf, idssb(j)+nres, iret)
 748             call xdrfint(ixdrf, jdssb(j)+nres, iret)
 749            else
 750             call xdrfint(ixdrf, ihpb(j), iret)
 751             call xdrfint(ixdrf, jhpb(j), iret)
 752            endif
 753       enddo
 754       call xdrffloat(ixdrf,real(eini),iret)
 755       call xdrffloat(ixdrf,real(efree),iret)
 756       call xdrffloat(ixdrf,real(rmsdev),iret)
 757       call xdrfint(ixdrf,iscor,iret)
 758 #endif
 759
 760       return
 761       end
 762 c------------------------------------------------------------------------------
 763       logical function conf_check(ii,iprint)
 764       implicit none
 765       include "DIMENSIONS"
 766       include "DIMENSIONS.ZSCOPT"
 767       include "DIMENSIONS.FREE"
 768 #ifdef MPI
 769       include "mpif.h"
 770       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
 771       include "COMMON.MPI"
 772 #endif
 773       include "COMMON.CHAIN"
 774       include "COMMON.IOUNITS"
 775       include "COMMON.PROTFILES"
 776       include "COMMON.NAMES"
 777       include "COMMON.VAR"
 778       include "COMMON.SBRIDGE"
 779       include "COMMON.GEO"
 780       include "COMMON.FFIELD"
 781       include "COMMON.ENEPS"
 782       include "COMMON.LOCAL"
 783       include "COMMON.WEIGHTS"
 784       include "COMMON.INTERACT"
 785       include "COMMON.FREE"
 786       include "COMMON.ENERGIES"
 787       include "COMMON.CONTROL"
 788       include "COMMON.TORCNSTR"
 789       integer j,k,l,ii,itj,iprint
 790 c      if (.not.check_conf) then
 791 c        conf_check=.true.
 792 c        return
 793 c      endif
 794       call int_from_cart1(.false.)
 795       do j=nnt+1,nct
 796         if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and.
 797      &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.5d0)) then
 798           if (iprint.gt.0)
 799      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
 800      &      " for conformation",ii
 801           if (iprint.gt.1) then
 802             write (iout,*) "The Cartesian geometry is:"
 803             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 804             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 805             write (iout,*) "The internal geometry is:"
 806             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 807             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 808             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 809             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 810             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 811             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 812           endif
 813           if (iprint.gt.0) write (iout,*)
 814      &      "This conformation WILL NOT be added to the database."
 815           conf_check=.false.
 816           return
 817         endif
 818       enddo
 819       do j=nnt,nct
 820         itj=itype(j)
 821         if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and.
 822      &     (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
 823           if (iprint.gt.0)
 824      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
 825      &     restyp(itj),itj,dsc(iabs(itj))," for conformation",ii
 826           if (iprint.gt.1) then
 827             write (iout,*) "The Cartesian geometry is:"
 828             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 829             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 830             write (iout,*) "The internal geometry is:"
 831             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 832             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 833             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 834             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 835             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 836             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 837           endif
 838           if (iprint.gt.0) write (iout,*)
 839      &      "This conformation WILL NOT be added to the database."
 840           conf_check=.false.
 841           return
 842         endif
 843       enddo
 844       do j=3,nres
 845         if (theta(j).le.0.0d0) then
 846           if (iprint.gt.0)
 847      &    write (iout,*) "Zero theta angle(s) in conformation",ii
 848           if (iprint.gt.1) then
 849             write (iout,*) "The Cartesian geometry is:"
 850             write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 851             write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 852             write (iout,*) "The internal geometry is:"
 853             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 854             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 855             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 856             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 857             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 858             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 859           endif
 860           if (iprint.gt.0) write (iout,*)
 861      &      "This conformation WILL NOT be added to the database."
 862           conf_check=.false.
 863           return
 864         endif
 865         if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
 866       enddo
 867       conf_check=.true.
 868 c      write (iout,*) "conf_check passed",ii
 869       return
 870       end