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