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