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