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