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