source/unres/src_MD/readrtns.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21       if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23       if (modecalc.eq.14) then
  24          call read_MDpar
  25          call read_REMDpar
  26       endif
  27 C Read MUCA control parameters if required
  28       if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31 csa      if (modecalc.eq.8) then
  32 csa       inquire (file="fort.40",exist=file_exist)
  33 csa       if (.not.file_exist) call csaread
  34 csa      endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39 C Print restraint information
  40 #ifdef MPI
  41       if (.not. out1file .or. me.eq.king) then
  42 #endif
  43       if (nhpb.gt.nss)
  44      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  45      & " distance constraints have been imposed"
  46       do i=nss+1,nhpb
  47         write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
  48      &     ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
  49       enddo
  50 #ifdef MPI
  51       endif
  52 #endif
  53 c      print *,"Processor",myrank," leaves READRTNS"
  54       return
  55       end
  56 C-------------------------------------------------------------------------------
  57       subroutine read_control
  58 C
  59 C Read contorl data
  60 C
  61       implicit real*8 (a-h,o-z)
  62       include 'DIMENSIONS'
  63 #ifdef MP
  64       include 'mpif.h'
  65       logical OKRandom, prng_restart
  66       real*8  r1
  67 #endif
  68       include 'COMMON.IOUNITS'
  69       include 'COMMON.TIME1'
  70       include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74       include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76 csa      include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78       include 'COMMON.MUCA'
  79       include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 100       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 101       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 102       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 103       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 104       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 105       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 106       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 107       call reada(controlcard,'DRMS',drms,0.1D0)
 108       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 109        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 110        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 111        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 112        write (iout,'(a,f10.1)')'DRMS    = ',drms
 113        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 114        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 115       endif
 116       call readi(controlcard,'NZ_START',nz_start,0)
 117       call readi(controlcard,'NZ_END',nz_end,0)
 118       call readi(controlcard,'IZ_SC',iz_sc,0)
 119       timlim=60.0D0*timlim
 120       safety = 60.0d0*safety
 121       timem=timlim
 122       modecalc=0
 123       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 124       minim=(index(controlcard,'MINIMIZE').gt.0)
 125       dccart=(index(controlcard,'CART').gt.0)
 126       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 127       overlapsc=.not.overlapsc
 128       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 129       searchsc=.not.searchsc
 130       sideadd=(index(controlcard,'SIDEADD').gt.0)
 131       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 132       outpdb=(index(controlcard,'PDBOUT').gt.0)
 133       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 134       pdbref=(index(controlcard,'PDBREF').gt.0)
 135       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 136       indpdb=index(controlcard,'PDBSTART')
 137       extconf=(index(controlcard,'EXTCONF').gt.0)
 138       call readi(controlcard,'IPRINT',iprint,0)
 139       call readi(controlcard,'MAXGEN',maxgen,10000)
 140       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 141       call readi(controlcard,"KDIAG",kdiag,0)
 142       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 143       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 144      & write (iout,*) "RESCALE_MODE",rescale_mode
 145       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 146       if (index(controlcard,'REGULAR').gt.0.0D0) then
 147         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 148         modecalc=1
 149         refstr=.true.
 150       endif
 151       if (index(controlcard,'CHECKGRAD').gt.0) then
 152         modecalc=5
 153         if (index(controlcard,'CART').gt.0) then
 154           icheckgrad=1
 155         elseif (index(controlcard,'CARINT').gt.0) then
 156           icheckgrad=2
 157         else
 158           icheckgrad=3
 159         endif
 160       elseif (index(controlcard,'THREAD').gt.0) then
 161         modecalc=2
 162         call readi(controlcard,'THREAD',nthread,0)
 163         if (nthread.gt.0) then
 164           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 165         else
 166           if (fg_rank.eq.0)
 167      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 168           stop 'Error termination in Read_Control.'
 169         endif
 170       else if (index(controlcard,'MCMA').gt.0) then
 171         modecalc=3
 172       else if (index(controlcard,'MCEE').gt.0) then
 173         modecalc=6
 174       else if (index(controlcard,'MULTCONF').gt.0) then
 175         modecalc=4
 176       else if (index(controlcard,'MAP').gt.0) then
 177         modecalc=7
 178         call readi(controlcard,'MAP',nmap,0)
 179       else if (index(controlcard,'CSA').gt.0) then
 180            write(*,*) "CSA not supported in this version"
 181            stop
 182 csa        modecalc=8
 183 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 184 crc
 185 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 186 crc
 187 crc        modecalc=9
 188 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 189 cfmc        modecalc=10
 190       else if (index(controlcard,'SOFTREG').gt.0) then
 191         modecalc=11
 192       else if (index(controlcard,'CHECK_BOND').gt.0) then
 193         modecalc=-1
 194       else if (index(controlcard,'TEST').gt.0) then
 195         modecalc=-2
 196       else if (index(controlcard,'MD').gt.0) then
 197         modecalc=12
 198       else if (index(controlcard,'RE ').gt.0) then
 199         modecalc=14
 200       endif
 201
 202       lmuca=index(controlcard,'MUCA').gt.0
 203       call readi(controlcard,'MUCADYN',mucadyn,0)
 204       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 205       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 206      & then
 207        write (iout,*) 'MUCADYN=',mucadyn
 208        write (iout,*) 'MUCASMOOTH=',muca_smooth
 209       endif
 210
 211       iscode=index(controlcard,'ONE_LETTER')
 212       indphi=index(controlcard,'PHI')
 213       indback=index(controlcard,'BACK')
 214       iranconf=index(controlcard,'RAND_CONF')
 215       i2ndstr=index(controlcard,'USE_SEC_PRED')
 216       gradout=index(controlcard,'GRADOUT').gt.0
 217       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 218
 219       if(me.eq.king.or..not.out1file)
 220      & write (iout,'(2a)') diagmeth(kdiag),
 221      &  ' routine used to diagonalize matrices.'
 222       return
 223       end
 224 c--------------------------------------------------------------------------
 225       subroutine read_REMDpar
 226 C
 227 C Read REMD settings
 228 C
 229       implicit real*8 (a-h,o-z)
 230       include 'DIMENSIONS'
 231       include 'COMMON.IOUNITS'
 232       include 'COMMON.TIME1'
 233       include 'COMMON.MD'
 234 #ifndef LANG0
 235       include 'COMMON.LANGEVIN'
 236 #else
 237       include 'COMMON.LANGEVIN.lang0'
 238 #endif
 239       include 'COMMON.INTERACT'
 240       include 'COMMON.NAMES'
 241       include 'COMMON.GEO'
 242       include 'COMMON.REMD'
 243       include 'COMMON.CONTROL'
 244       include 'COMMON.SETUP'
 245       character*80 ucase
 246       character*320 controlcard
 247       character*3200 controlcard1
 248       integer iremd_m_total
 249
 250       if(me.eq.king.or..not.out1file)
 251      & write (iout,*) "REMD setup"
 252
 253       call card_concat(controlcard)
 254       call readi(controlcard,"NREP",nrep,3)
 255       call readi(controlcard,"NSTEX",nstex,1000)
 256       call reada(controlcard,"RETMIN",retmin,10.0d0)
 257       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 258       mremdsync=(index(controlcard,'SYNC').gt.0)
 259       call readi(controlcard,"NSYN",i_sync_step,100)
 260       restart1file=(index(controlcard,'REST1FILE').gt.0)
 261       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 262       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 263       if(max_cache_traj_use.gt.max_cache_traj)
 264      &           max_cache_traj_use=max_cache_traj
 265       if(me.eq.king.or..not.out1file) then
 266 cd       if (traj1file) then
 267 crc caching is in testing - NTWX is not ignored
 268 cd        write (iout,*) "NTWX value is ignored"
 269 cd        write (iout,*) "  trajectory is stored to one file by master"
 270 cd        write (iout,*) "  before exchange at NSTEX intervals"
 271 cd       endif
 272        write (iout,*) "NREP= ",nrep
 273        write (iout,*) "NSTEX= ",nstex
 274        write (iout,*) "SYNC= ",mremdsync
 275        write (iout,*) "NSYN= ",i_sync_step
 276        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 277       endif
 278
 279       t_exchange_only=(index(controlcard,'TONLY').gt.0)
 280       call readi(controlcard,"HREMD",hremd,0)
 281       if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
 282         write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
 283       endif
 284       if(usampl.and.hremd.gt.0) then
 285             write (iout,'(//a)')
 286      &      "========== ERROR: USAMPL and HREMD cannot be used together"
 287 #ifdef MPI
 288             call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 289 #endif
 290             stop
 291       endif
 292
 293
 294       remd_tlist=.false.
 295       if (index(controlcard,'TLIST').gt.0) then
 296          remd_tlist=.true.
 297          call card_concat(controlcard1)
 298          read(controlcard1,*) (remd_t(i),i=1,nrep)
 299          if(me.eq.king.or..not.out1file)
 300      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 301       endif
 302       remd_mlist=.false.
 303       if (index(controlcard,'MLIST').gt.0) then
 304          remd_mlist=.true.
 305          call card_concat(controlcard1)
 306          read(controlcard1,*) (remd_m(i),i=1,nrep)
 307          if(me.eq.king.or..not.out1file) then
 308           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 309           iremd_m_total=0
 310           do i=1,nrep
 311            iremd_m_total=iremd_m_total+remd_m(i)
 312           enddo
 313           if(hremd.gt.1)then
 314            write (iout,*) 'Total number of replicas ',
 315      &       iremd_m_total*hremd
 316           else
 317            write (iout,*) 'Total number of replicas ',iremd_m_total
 318           endif
 319          endif
 320       endif
 321       if(me.eq.king.or..not.out1file)
 322      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 323       return
 324       end
 325 c--------------------------------------------------------------------------
 326       subroutine read_MDpar
 327 C
 328 C Read MD settings
 329 C
 330       implicit real*8 (a-h,o-z)
 331       include 'DIMENSIONS'
 332       include 'COMMON.IOUNITS'
 333       include 'COMMON.TIME1'
 334       include 'COMMON.MD'
 335 #ifndef LANG0
 336       include 'COMMON.LANGEVIN'
 337 #else
 338       include 'COMMON.LANGEVIN.lang0'
 339 #endif
 340       include 'COMMON.INTERACT'
 341       include 'COMMON.NAMES'
 342       include 'COMMON.GEO'
 343       include 'COMMON.SETUP'
 344       include 'COMMON.CONTROL'
 345       include 'COMMON.SPLITELE'
 346       character*80 ucase
 347       character*320 controlcard
 348
 349       call card_concat(controlcard)
 350       call readi(controlcard,"NSTEP",n_timestep,1000000)
 351       call readi(controlcard,"NTWE",ntwe,100)
 352       call readi(controlcard,"NTWX",ntwx,1000)
 353       call reada(controlcard,"DT",d_time,1.0d-1)
 354       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 355       call reada(controlcard,"DAMAX",damax,1.0d1)
 356       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 357       call readi(controlcard,"LANG",lang,0)
 358       RESPA = index(controlcard,"RESPA") .gt. 0
 359       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 360       ntime_split0=ntime_split
 361       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 362       ntime_split0=ntime_split
 363       call reada(controlcard,"R_CUT",r_cut,2.0d0)
 364       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 365       rest = index(controlcard,"REST").gt.0
 366       tbf = index(controlcard,"TBF").gt.0
 367       call readi(controlcard,"HMC",hmc,0)
 368       tnp = index(controlcard,"NOSEPOINCARE99").gt.0
 369       tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
 370       tnh = index(controlcard,"NOSEHOOVER96").gt.0
 371       if (RESPA.and.tnh)then
 372         xiresp = index(controlcard,"XIRESP").gt.0
 373       endif
 374       call reada(controlcard,"Q_NP",Q_np,0.1d0)
 375       usampl = index(controlcard,"USAMPL").gt.0
 376
 377       mdpdb = index(controlcard,"MDPDB").gt.0
 378       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 379       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 380       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 381       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 382       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 383       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 384       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 385       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 386       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 387       large = index(controlcard,"LARGE").gt.0
 388       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 389       rattle = index(controlcard,"RATTLE").gt.0
 390 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 391       nset=0
 392       if(usampl) then
 393         call read_fragments
 394       endif
 395
 396       if(me.eq.king.or..not.out1file) then
 397        write (iout,*)
 398        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 399        write (iout,*)
 400        write (iout,'(a)') "The units are:"
 401        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 402        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 403      &  " acceleration: angstrom/(48.9 fs)**2"
 404        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 405        write (iout,*)
 406        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 407        write (iout,'(a60,f10.5,a)')
 408      &  "Initial time step of numerical integration:",d_time,
 409      &  " natural units"
 410        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 411        if (RESPA) then
 412         write (iout,'(2a,i4,a)')
 413      &    "A-MTS algorithm used; initial time step for fast-varying",
 414      &    " short-range forces split into",ntime_split," steps."
 415         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 416      &   r_cut," lambda",rlamb
 417        endif
 418        write (iout,'(2a,f10.5)')
 419      &  "Maximum acceleration threshold to reduce the time step",
 420      &  "/increase split number:",damax
 421        write (iout,'(2a,f10.5)')
 422      &  "Maximum predicted energy drift to reduce the timestep",
 423      &  "/increase split number:",edriftmax
 424        write (iout,'(a60,f10.5)')
 425      & "Maximum velocity threshold to reduce velocities:",dvmax
 426        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 427        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 428        if (rattle) write (iout,'(a60)')
 429      &  "Rattle algorithm used to constrain the virtual bonds"
 430       endif
 431       reset_fricmat=1000
 432       if (lang.gt.0) then
 433         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 434         call reada(controlcard,"RWAT",rwat,1.4d0)
 435         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 436         surfarea=index(controlcard,"SURFAREA").gt.0
 437         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 438         if(me.eq.king.or..not.out1file)then
 439          write (iout,'(/a,$)') "Langevin dynamics calculation"
 440          if (lang.eq.1) then
 441           write (iout,'(a/)')
 442      &      " with direct integration of Langevin equations"
 443          else if (lang.eq.2) then
 444           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 445          else if (lang.eq.3) then
 446           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 447          else if (lang.eq.4) then
 448           write (iout,'(a/)') " in overdamped mode"
 449          else
 450           write (iout,'(//a,i5)')
 451      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 452           stop
 453          endif
 454          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 455          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 456          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 457          write (iout,'(a60,f10.5)')
 458      &   "Scaling factor of the friction forces:",scal_fric
 459          if (surfarea) write (iout,'(2a,i10,a)')
 460      &     "Friction coefficients will be scaled by solvent-accessible",
 461      &     " surface area every",reset_fricmat," steps."
 462         endif
 463 c Calculate friction coefficients and bounds of stochastic forces
 464         eta=6*pi*cPoise*etawat
 465         if(me.eq.king.or..not.out1file)
 466      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 467      &   ,eta
 468         gamp=scal_fric*(pstok+rwat)*eta
 469         stdfp=dsqrt(2*Rb*t_bath/d_time)
 470         do i=1,ntyp
 471           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 472           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 473         enddo
 474         if(me.eq.king.or..not.out1file)then
 475          write (iout,'(/2a/)')
 476      &   "Radii of site types and friction coefficients and std's of",
 477      &   " stochastic forces of fully exposed sites"
 478          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 479          do i=1,ntyp
 480           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 481      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 482          enddo
 483         endif
 484       else if (tbf) then
 485         if(me.eq.king.or..not.out1file)then
 486          write (iout,'(a)') "Berendsen bath calculation"
 487          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 488          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 489          if (reset_moment)
 490      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 491      &   count_reset_moment," steps"
 492          if (reset_vel)
 493      &    write (iout,'(a,i10,a)')
 494      &    "Velocities will be reset at random every",count_reset_vel,
 495      &   " steps"
 496         endif
 497       else if (tnp .or. tnp1 .or. tnh) then
 498         if (tnp .or. tnp1) then
 499            write (iout,'(a)') "Nose-Poincare bath calculation"
 500            if (tnp) write (iout,'(a)')
 501      & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
 502            if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
 503         else
 504            write (iout,'(a)') "Nose-Hoover bath calculation"
 505            write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
 506               nresn=1
 507               nyosh=1
 508               nnos=1
 509               do i=1,nnos
 510                qmass(i)=Q_np
 511                xlogs(i)=1.0
 512                vlogs(i)=0.0
 513               enddo
 514               do i=1,nyosh
 515                WDTI(i) = 1.0*d_time/nresn
 516                WDTI2(i)=WDTI(i)/2
 517                WDTI4(i)=WDTI(i)/4
 518                WDTI8(i)=WDTI(i)/8
 519               enddo
 520               if (RESPA) then
 521                if(xiresp) then
 522                  write (iout,'(a)') "NVT-XI-RESPA algorithm"
 523                else
 524                  write (iout,'(a)') "NVT-XO-RESPA algorithm"
 525                endif
 526                do i=1,nyosh
 527                 WDTIi(i) = 1.0*d_time/nresn/ntime_split
 528                 WDTIi2(i)=WDTIi(i)/2
 529                 WDTIi4(i)=WDTIi(i)/4
 530                 WDTIi8(i)=WDTIi(i)/8
 531                enddo
 532               endif
 533         endif
 534
 535         write (iout,'(a60,f10.5)') "Temperature:",t_bath
 536         write (iout,'(a60,f10.5)') "Q =",Q_np
 537         if (reset_moment)
 538      &  write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 539      &   count_reset_moment," steps"
 540         if (reset_vel)
 541      &    write (iout,'(a,i10,a)')
 542      &    "Velocities will be reset at random every",count_reset_vel,
 543      &   " steps"
 544
 545       else if (hmc.gt.0) then
 546          write (iout,'(a)') "Hybrid Monte Carlo calculation"
 547          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 548          write (iout,'(a60,i10)')
 549      &         "Number of MD steps between Metropolis tests:",hmc
 550
 551       else
 552         if(me.eq.king.or..not.out1file)
 553      &   write (iout,'(a31)') "Microcanonical mode calculation"
 554       endif
 555       if(me.eq.king.or..not.out1file)then
 556        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 557        if (usampl) then
 558           write(iout,*) "MD running with constraints."
 559           write(iout,*) "Equilibration time ", eq_time, " mtus."
 560           write(iout,*) "Constraining ", nfrag," fragments."
 561           write(iout,*) "Length of each fragment, weight and q0:"
 562           do iset=1,nset
 563            write (iout,*) "Set of restraints #",iset
 564            do i=1,nfrag
 565               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 566      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 567            enddo
 568            write(iout,*) "constraints between ", npair, "fragments."
 569            write(iout,*) "constraint pairs, weights and q0:"
 570            do i=1,npair
 571             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 572      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 573            enddo
 574            write(iout,*) "angle constraints within ", nfrag_back,
 575      &      "backbone fragments."
 576            write(iout,*) "fragment, weights:"
 577            do i=1,nfrag_back
 578             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 579      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 580      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 581            enddo
 582           enddo
 583         iset=mod(kolor,nset)+1
 584        endif
 585       endif
 586       if(me.eq.king.or..not.out1file)
 587      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 588       return
 589       end
 590 c------------------------------------------------------------------------------
 591       subroutine molread
 592 C
 593 C Read molecular data.
 594 C
 595       implicit real*8 (a-h,o-z)
 596       include 'DIMENSIONS'
 597 #ifdef MPI
 598       include 'mpif.h'
 599       integer error_msg
 600 #endif
 601       include 'COMMON.IOUNITS'
 602       include 'COMMON.GEO'
 603       include 'COMMON.VAR'
 604       include 'COMMON.INTERACT'
 605       include 'COMMON.LOCAL'
 606       include 'COMMON.NAMES'
 607       include 'COMMON.CHAIN'
 608       include 'COMMON.FFIELD'
 609       include 'COMMON.SBRIDGE'
 610       include 'COMMON.HEADER'
 611       include 'COMMON.CONTROL'
 612       include 'COMMON.DBASE'
 613       include 'COMMON.THREAD'
 614       include 'COMMON.CONTACTS'
 615       include 'COMMON.TORCNSTR'
 616       include 'COMMON.TIME1'
 617       include 'COMMON.BOUNDS'
 618       include 'COMMON.MD'
 619       include 'COMMON.REMD'
 620       include 'COMMON.SETUP'
 621       character*4 sequence(maxres)
 622       integer rescode
 623       double precision x(maxvar)
 624       character*256 pdbfile
 625       character*320 weightcard
 626       character*80 weightcard_t,ucase
 627       dimension itype_pdb(maxres)
 628       common /pizda/ itype_pdb
 629       logical seq_comp,fail
 630       double precision energia(0:n_ene)
 631       integer ilen
 632       external ilen
 633 C
 634 C Body
 635 C
 636 C Read weights of the subsequent energy terms.
 637       if(hremd.gt.0) then
 638
 639        k=0
 640        do il=1,hremd
 641         do i=1,nrep
 642          do j=1,remd_m(i)
 643           i2set(k)=il
 644           k=k+1
 645          enddo
 646         enddo
 647        enddo
 648
 649        if(me.eq.king.or..not.out1file) then
 650         write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
 651         write (iout,*) 'Current weights for processor ',
 652      &                 me,' set ',i2set(me)
 653        endif
 654
 655        do i=1,hremd
 656          call card_concat(weightcard)
 657          call reada(weightcard,'WLONG',wlong,1.0D0)
 658          call reada(weightcard,'WSC',wsc,wlong)
 659          call reada(weightcard,'WSCP',wscp,wlong)
 660          call reada(weightcard,'WELEC',welec,1.0D0)
 661          call reada(weightcard,'WVDWPP',wvdwpp,welec)
 662          call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 663          call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 664          call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 665          call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 666          call reada(weightcard,'WTURN3',wturn3,1.0D0)
 667          call reada(weightcard,'WTURN4',wturn4,1.0D0)
 668          call reada(weightcard,'WTURN6',wturn6,1.0D0)
 669          call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 670          call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 671          call reada(weightcard,'WBOND',wbond,1.0D0)
 672          call reada(weightcard,'WTOR',wtor,1.0D0)
 673          call reada(weightcard,'WTORD',wtor_d,1.0D0)
 674          call reada(weightcard,'WANG',wang,1.0D0)
 675          call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 676          call reada(weightcard,'SCAL14',scal14,0.4D0)
 677          call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 678          call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 679          call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 680          call reada(weightcard,'TEMP0',temp0,300.0d0)
 681          if (index(weightcard,'SOFT').gt.0) ipot=6
 682 C 12/1/95 Added weight for the multi-body term WCORR
 683          call reada(weightcard,'WCORRH',wcorr,1.0D0)
 684          if (wcorr4.gt.0.0d0) wcorr=wcorr4
 685
 686          hweights(i,1)=wsc
 687          hweights(i,2)=wscp
 688          hweights(i,3)=welec
 689          hweights(i,4)=wcorr
 690          hweights(i,5)=wcorr5
 691          hweights(i,6)=wcorr6
 692          hweights(i,7)=wel_loc
 693          hweights(i,8)=wturn3
 694          hweights(i,9)=wturn4
 695          hweights(i,10)=wturn6
 696          hweights(i,11)=wang
 697          hweights(i,12)=wscloc
 698          hweights(i,13)=wtor
 699          hweights(i,14)=wtor_d
 700          hweights(i,15)=wstrain
 701          hweights(i,16)=wvdwpp
 702          hweights(i,17)=wbond
 703          hweights(i,18)=scal14
 704          hweights(i,21)=wsccor
 705
 706        enddo
 707
 708        do i=1,n_ene
 709          weights(i)=hweights(i2set(me),i)
 710        enddo
 711        wsc    =weights(1)
 712        wscp   =weights(2)
 713        welec  =weights(3)
 714        wcorr  =weights(4)
 715        wcorr5 =weights(5)
 716        wcorr6 =weights(6)
 717        wel_loc=weights(7)
 718        wturn3 =weights(8)
 719        wturn4 =weights(9)
 720        wturn6 =weights(10)
 721        wang   =weights(11)
 722        wscloc =weights(12)
 723        wtor   =weights(13)
 724        wtor_d =weights(14)
 725        wstrain=weights(15)
 726        wvdwpp =weights(16)
 727        wbond  =weights(17)
 728        scal14 =weights(18)
 729        wsccor =weights(21)
 730
 731
 732       else
 733        call card_concat(weightcard)
 734        call reada(weightcard,'WLONG',wlong,1.0D0)
 735        call reada(weightcard,'WSC',wsc,wlong)
 736        call reada(weightcard,'WSCP',wscp,wlong)
 737        call reada(weightcard,'WELEC',welec,1.0D0)
 738        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 739        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 740        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 741        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 742        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 743        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 744        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 745        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 746        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 747        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 748        call reada(weightcard,'WBOND',wbond,1.0D0)
 749        call reada(weightcard,'WTOR',wtor,1.0D0)
 750        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 751        call reada(weightcard,'WANG',wang,1.0D0)
 752        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 753        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 754        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 755        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 756        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 757        call reada(weightcard,'SCAL14',scal14,0.4D0)
 758        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 759        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 760        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 761        call reada(weightcard,'TEMP0',temp0,300.0d0)
 762        if (index(weightcard,'SOFT').gt.0) ipot=6
 763 C 12/1/95 Added weight for the multi-body term WCORR
 764        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 765        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 766        weights(1)=wsc
 767        weights(2)=wscp
 768        weights(3)=welec
 769        weights(4)=wcorr
 770        weights(5)=wcorr5
 771        weights(6)=wcorr6
 772        weights(7)=wel_loc
 773        weights(8)=wturn3
 774        weights(9)=wturn4
 775        weights(10)=wturn6
 776        weights(11)=wang
 777        weights(12)=wscloc
 778        weights(13)=wtor
 779        weights(14)=wtor_d
 780        weights(15)=wstrain
 781        weights(16)=wvdwpp
 782        weights(17)=wbond
 783        weights(18)=scal14
 784        weights(21)=wsccor
 785       endif
 786        weights(25)=wdfa_dist
 787        weights(26)=wdfa_tor
 788        weights(27)=wdfa_nei
 789        weights(28)=wdfa_beta
 790
 791       if(me.eq.king.or..not.out1file)
 792      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 793      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 794      &  wturn4,wturn6,
 795      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 796    10 format (/'Energy-term weights (unscaled):'//
 797      & 'WSCC=   ',f10.6,' (SC-SC)'/
 798      & 'WSCP=   ',f10.6,' (SC-p)'/
 799      & 'WELEC=  ',f10.6,' (p-p electr)'/
 800      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 801      & 'WBOND=  ',f10.6,' (stretching)'/
 802      & 'WANG=   ',f10.6,' (bending)'/
 803      & 'WSCLOC= ',f10.6,' (SC local)'/
 804      & 'WTOR=   ',f10.6,' (torsional)'/
 805      & 'WTORD=  ',f10.6,' (double torsional)'/
 806      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 807      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 808      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 809      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 810      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 811      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 812      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 813      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 814      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 815      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 816      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 817      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 818      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 819       if(me.eq.king.or..not.out1file)then
 820        if (wcorr4.gt.0.0d0) then
 821         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 822      &   'between contact pairs of peptide groups'
 823         write (iout,'(2(a,f5.3/))')
 824      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 825      &  'Range of quenching the correlation terms:',2*delt_corr
 826        else if (wcorr.gt.0.0d0) then
 827         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 828      &   'between contact pairs of peptide groups'
 829        endif
 830        write (iout,'(a,f8.3)')
 831      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 832        write (iout,'(a,f8.3)')
 833      &  'General scaling factor of SC-p interactions:',scalscp
 834       endif
 835       r0_corr=cutoff_corr-delt_corr
 836       do i=1,20
 837         aad(i,1)=scalscp*aad(i,1)
 838         aad(i,2)=scalscp*aad(i,2)
 839         bad(i,1)=scalscp*bad(i,1)
 840         bad(i,2)=scalscp*bad(i,2)
 841       enddo
 842       call rescale_weights(t_bath)
 843       if(me.eq.king.or..not.out1file)
 844      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 845      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 846      &  wturn4,wturn6,
 847      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 848    22 format (/'Energy-term weights (scaled):'//
 849      & 'WSCC=   ',f10.6,' (SC-SC)'/
 850      & 'WSCP=   ',f10.6,' (SC-p)'/
 851      & 'WELEC=  ',f10.6,' (p-p electr)'/
 852      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 853      & 'WBOND=  ',f10.6,' (stretching)'/
 854      & 'WANG=   ',f10.6,' (bending)'/
 855      & 'WSCLOC= ',f10.6,' (SC local)'/
 856      & 'WTOR=   ',f10.6,' (torsional)'/
 857      & 'WTORD=  ',f10.6,' (double torsional)'/
 858      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 859      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 860      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 861      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 862      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 863      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 864      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 865      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 866      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 867      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 868      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 869      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 870      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 871       if(me.eq.king.or..not.out1file)
 872      & write (iout,*) "Reference temperature for weights calculation:",
 873      &  temp0
 874       call reada(weightcard,"D0CM",d0cm,3.78d0)
 875       call reada(weightcard,"AKCM",akcm,15.1d0)
 876       call reada(weightcard,"AKTH",akth,11.0d0)
 877       call reada(weightcard,"AKCT",akct,12.0d0)
 878       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 879       call reada(weightcard,"V2SS",v2ss,7.61d0)
 880       call reada(weightcard,"V3SS",v3ss,13.7d0)
 881       call reada(weightcard,"EBR",ebr,-5.50D0)
 882       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 883       do i=1,maxres
 884         dyn_ss_mask(i)=.false.
 885       enddo
 886       do i=1,maxres-1
 887         do j=i+1,maxres
 888           dyn_ssbond_ij(i,j)=1.0d300
 889         enddo
 890       enddo
 891       call reada(weightcard,"HT",Ht,0.0D0)
 892       if (dyn_ss) then
 893         ss_depth=ebr/wsc-0.25*eps(1,1)
 894         Ht=Ht/wsc-0.25*eps(1,1)
 895         akcm=akcm*wstrain/wsc
 896         akth=akth*wstrain/wsc
 897         akct=akct*wstrain/wsc
 898         v1ss=v1ss*wstrain/wsc
 899         v2ss=v2ss*wstrain/wsc
 900         v3ss=v3ss*wstrain/wsc
 901       else
 902         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 903       endif
 904
 905       if(me.eq.king.or..not.out1file) then
 906        write (iout,*) "Parameters of the SS-bond potential:"
 907        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 908      & " AKCT",akct
 909        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 910        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 911        write (iout,*)" HT",Ht
 912        print *,'indpdb=',indpdb,' pdbref=',pdbref
 913       endif
 914       if (indpdb.gt.0 .or. pdbref) then
 915         read(inp,'(a)') pdbfile
 916         if(me.eq.king.or..not.out1file)
 917      &   write (iout,'(2a)') 'PDB data will be read from file ',
 918      &   pdbfile(:ilen(pdbfile))
 919         open(ipdbin,file=pdbfile,status='old',err=33)
 920         goto 34
 921   33    write (iout,'(a)') 'Error opening PDB file.'
 922         stop
 923   34    continue
 924 c        print *,'Begin reading pdb data'
 925         call readpdb
 926 c        print *,'Finished reading pdb data'
 927         if(me.eq.king.or..not.out1file)
 928      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 929      &   ' nstart_sup=',nstart_sup
 930         do i=1,nres
 931           itype_pdb(i)=itype(i)
 932         enddo
 933         close (ipdbin)
 934         nnt=nstart_sup
 935         nct=nstart_sup+nsup-1
 936         call contact(.false.,ncont_ref,icont_ref,co)
 937
 938         if (sideadd) then
 939 C Following 2 lines for diagnostics; comment out if not needed
 940          write (iout,*) "Before sideadd"
 941          call intout
 942          if(me.eq.king.or..not.out1file)
 943      &    write(iout,*)'Adding sidechains'
 944          maxsi=1000
 945          do i=2,nres-1
 946           iti=itype(i)
 947           if (iti.ne.10) then
 948             nsi=0
 949             fail=.true.
 950             do while (fail.and.nsi.le.maxsi)
 951               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 952               nsi=nsi+1
 953             enddo
 954             if(fail) write(iout,*)'Adding sidechain failed for res ',
 955      &              i,' after ',nsi,' trials'
 956           endif
 957          enddo
 958 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
 959          call chainbuild
 960         endif
 961 C Following 2 lines for diagnostics; comment out if not needed
 962 c        write (iout,*) "After sideadd"
 963 c        call intout
 964       endif
 965       if (indpdb.eq.0) then
 966 C Read sequence if not taken from the pdb file.
 967         read (inp,*) nres
 968 c        print *,'nres=',nres
 969         if (iscode.gt.0) then
 970           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 971         else
 972           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 973         endif
 974 C Convert sequence to numeric code
 975         do i=1,nres
 976           itype(i)=rescode(i,sequence(i),iscode)
 977         enddo
 978 C Assign initial virtual bond lengths
 979         do i=2,nres
 980           vbld(i)=vbl
 981           vbld_inv(i)=vblinv
 982         enddo
 983         do i=2,nres-1
 984           vbld(i+nres)=dsc(itype(i))
 985           vbld_inv(i+nres)=dsc_inv(itype(i))
 986 c          write (iout,*) "i",i," itype",itype(i),
 987 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 988         enddo
 989       endif
 990 c      print *,nres
 991 c      print '(20i4)',(itype(i),i=1,nres)
 992       do i=1,nres
 993 #ifdef PROCOR
 994         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 995 #else
 996         if (itype(i).eq.21) then
 997 #endif
 998           itel(i)=0
 999 #ifdef PROCOR
1000         else if (itype(i+1).ne.20) then
1001 #else
1002         else if (itype(i).ne.20) then
1003 #endif
1004           itel(i)=1
1005         else
1006           itel(i)=2
1007         endif
1008       enddo
1009       if(me.eq.king.or..not.out1file)then
1010        write (iout,*) "ITEL"
1011        do i=1,nres-1
1012          write (iout,*) i,itype(i),itel(i)
1013        enddo
1014        print *,'Call Read_Bridge.'
1015       endif
1016       call read_bridge
1017 C 8/13/98 Set limits to generating the dihedral angles
1018       do i=1,nres
1019         phibound(1,i)=-pi
1020         phibound(2,i)=pi
1021       enddo
1022       read (inp,*) ndih_constr
1023       if (ndih_constr.gt.0) then
1024         read (inp,*) ftors
1025         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1026         if(me.eq.king.or..not.out1file)then
1027          write (iout,*)
1028      &   'There are',ndih_constr,' constraints on phi angles.'
1029          do i=1,ndih_constr
1030           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1031          enddo
1032         endif
1033         do i=1,ndih_constr
1034           phi0(i)=deg2rad*phi0(i)
1035           drange(i)=deg2rad*drange(i)
1036         enddo
1037         if(me.eq.king.or..not.out1file)
1038      &   write (iout,*) 'FTORS',ftors
1039         do i=1,ndih_constr
1040           ii = idih_constr(i)
1041           phibound(1,ii) = phi0(i)-drange(i)
1042           phibound(2,ii) = phi0(i)+drange(i)
1043         enddo
1044       endif
1045       nnt=1
1046 #ifdef MPI
1047       if (me.eq.king) then
1048 #endif
1049        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1050        do i=1,nres
1051          write (iout,'(a3,i5,2f10.1)')
1052      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1053        enddo
1054 #ifdef MP
1055       endif
1056 #endif
1057       nct=nres
1058 cd      print *,'NNT=',NNT,' NCT=',NCT
1059       if (itype(1).eq.21) nnt=2
1060       if (itype(nres).eq.21) nct=nct-1
1061
1062 C     Bartek:READ init_vars
1063 C     Initialize variables!
1064 C     Juyong:READ read_info
1065 C     READ fragment information!!
1066 C     both routines should be in dfa.F file!!
1067
1068       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1069      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1070        call init_dfa_vars
1071        print*, 'init_dfa_vars finished!'
1072        call read_dfa_info
1073        print*, 'read_dfa_info finished!'
1074       endif
1075 C
1076       if (pdbref) then
1077         if(me.eq.king.or..not.out1file)
1078      &   write (iout,'(a,i3)') 'nsup=',nsup
1079         nstart_seq=nnt
1080         if (nsup.le.(nct-nnt+1)) then
1081           do i=0,nct-nnt+1-nsup
1082             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1083               nstart_seq=nnt+i
1084               goto 111
1085             endif
1086           enddo
1087           write (iout,'(a)')
1088      &            'Error - sequences to be superposed do not match.'
1089           stop
1090         else
1091           do i=0,nsup-(nct-nnt+1)
1092             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1093      &      then
1094               nstart_sup=nstart_sup+i
1095               nsup=nct-nnt+1
1096               goto 111
1097             endif
1098           enddo
1099           write (iout,'(a)')
1100      &            'Error - sequences to be superposed do not match.'
1101         endif
1102   111   continue
1103         if (nsup.eq.0) nsup=nct-nnt
1104         if (nstart_sup.eq.0) nstart_sup=nnt
1105         if (nstart_seq.eq.0) nstart_seq=nnt
1106         if(me.eq.king.or..not.out1file)
1107      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1108      &                 ' nstart_seq=',nstart_seq
1109       endif
1110 c--- Zscore rms -------
1111       if (nz_start.eq.0) nz_start=nnt
1112       if (nz_end.eq.0 .and. nsup.gt.0) then
1113         nz_end=nnt+nsup-1
1114       else if (nz_end.eq.0) then
1115         nz_end=nct
1116       endif
1117       if(me.eq.king.or..not.out1file)then
1118        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1119        write (iout,*) 'IZ_SC=',iz_sc
1120       endif
1121 c----------------------
1122       call init_int_table
1123       if (refstr) then
1124         if (.not.pdbref) then
1125           call read_angles(inp,*38)
1126           goto 39
1127    38     write (iout,'(a)') 'Error reading reference structure.'
1128 #ifdef MPI
1129           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1130           stop 'Error reading reference structure'
1131 #endif
1132    39     call chainbuild
1133           call setup_var
1134 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1135           nstart_sup=nnt
1136           nstart_seq=nnt
1137           nsup=nct-nnt+1
1138           do i=1,2*nres
1139             do j=1,3
1140               cref(j,i)=c(j,i)
1141             enddo
1142           enddo
1143           call contact(.true.,ncont_ref,icont_ref,co)
1144         endif
1145         if(me.eq.king.or..not.out1file)
1146      &   write (iout,*) 'Contact order:',co
1147         if (pdbref) then
1148         if(me.eq.king.or..not.out1file)
1149      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1150         do i=1,ncont_ref
1151           do j=1,2
1152             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1153           enddo
1154           if(me.eq.king.or..not.out1file)
1155      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1156      &     icont_ref(1,i),' ',
1157      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1158         enddo
1159         endif
1160       endif
1161 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1162       if (constr_dist.gt.0) then
1163         call read_dist_constr
1164       endif
1165
1166
1167       if (constr_homology.gt.0) then
1168         call read_constr_homology
1169       else
1170         homol_nset=0
1171       endif
1172
1173
1174       if (nhpb.gt.0) call hpb_partition
1175 c      write (iout,*) "After read_dist_constr nhpb",nhpb
1176 c      call flush(iout)
1177       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1178      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1179      &    modecalc.ne.10) then
1180 C If input structure hasn't been supplied from the PDB file read or generate
1181 C initial geometry.
1182         if (iranconf.eq.0 .and. .not. extconf) then
1183           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1184      &     write (iout,'(a)') 'Initial geometry will be read in.'
1185           if (read_cart) then
1186             read(inp,'(8f10.5)',end=36,err=36)
1187      &       ((c(l,k),l=1,3),k=1,nres),
1188      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1189             call int_from_cart1(.false.)
1190             do i=1,nres-1
1191               do j=1,3
1192                 dc(j,i)=c(j,i+1)-c(j,i)
1193                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1194               enddo
1195             enddo
1196             do i=nnt,nct
1197               if (itype(i).ne.10) then
1198                 do j=1,3
1199                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1200                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1201                 enddo
1202               endif
1203             enddo
1204             return
1205           else
1206             call read_angles(inp,*36)
1207           endif
1208           goto 37
1209    36     write (iout,'(a)') 'Error reading angle file.'
1210 #ifdef MPI
1211           call mpi_finalize( MPI_COMM_WORLD,IERR )
1212 #endif
1213           stop 'Error reading angle file.'
1214    37     continue
1215         else if (extconf) then
1216          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1217      &    write (iout,'(a)') 'Extended chain initial geometry.'
1218          do i=3,nres
1219           theta(i)=90d0*deg2rad
1220          enddo
1221          do i=4,nres
1222           phi(i)=180d0*deg2rad
1223          enddo
1224          do i=2,nres-1
1225           alph(i)=110d0*deg2rad
1226          enddo
1227          do i=2,nres-1
1228           omeg(i)=-120d0*deg2rad
1229          enddo
1230         else
1231           if(me.eq.king.or..not.out1file)
1232      &     write (iout,'(a)') 'Random-generated initial geometry.'
1233
1234
1235 #ifdef MPI
1236           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1237      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1238 #endif
1239             do itrial=1,100
1240               itmp=1
1241               call gen_rand_conf(itmp,*30)
1242               goto 40
1243    30         write (iout,*) 'Failed to generate random conformation',
1244      &          ', itrial=',itrial
1245               write (*,*) 'Processor:',me,
1246      &          ' Failed to generate random conformation',
1247      &          ' itrial=',itrial
1248               call intout
1249
1250 #ifdef AIX
1251               call flush_(iout)
1252 #else
1253               call flush(iout)
1254 #endif
1255             enddo
1256             write (iout,'(a,i3,a)') 'Processor:',me,
1257      &        ' error in generating random conformation.'
1258             write (*,'(a,i3,a)') 'Processor:',me,
1259      &        ' error in generating random conformation.'
1260             call flush(iout)
1261 #ifdef MPI
1262             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1263    40       continue
1264           endif
1265 #else
1266    40     continue
1267 #endif
1268         endif
1269       elseif (modecalc.eq.4) then
1270         read (inp,'(a)') intinname
1271         open (intin,file=intinname,status='old',err=333)
1272         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1273      &  write (iout,'(a)') 'intinname',intinname
1274         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1275         goto 334
1276   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1277 #ifdef MPI
1278         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1279 #endif
1280         stop 'Error opening angle file.'
1281   334   continue
1282
1283       endif
1284 C Generate distance constraints, if the PDB structure is to be regularized.
1285       if (nthread.gt.0) then
1286         call read_threadbase
1287       endif
1288       call setup_var
1289       if (me.eq.king .or. .not. out1file)
1290      & call intout
1291       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1292         write (iout,'(/a,i3,a)')
1293      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1294         write (iout,'(20i4)') (iss(i),i=1,ns)
1295        if (dyn_ss) then
1296           write(iout,*)"Running with dynamic disulfide-bond formation"
1297        else
1298         write (iout,'(/a/)') 'Pre-formed links are:'
1299         do i=1,nss
1300           i1=ihpb(i)-nres
1301           i2=jhpb(i)-nres
1302           it1=itype(i1)
1303           it2=itype(i2)
1304           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1305      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1306      &    ebr,forcon(i)
1307         enddo
1308         write (iout,'(a)')
1309        endif
1310       endif
1311       if (ns.gt.0.and.dyn_ss) then
1312           do i=nss+1,nhpb
1313             ihpb(i-nss)=ihpb(i)
1314             jhpb(i-nss)=jhpb(i)
1315             forcon(i-nss)=forcon(i)
1316             dhpb(i-nss)=dhpb(i)
1317           enddo
1318           nhpb=nhpb-nss
1319           nss=0
1320           call hpb_partition
1321           do i=1,ns
1322             dyn_ss_mask(iss(i))=.true.
1323           enddo
1324       endif
1325       if (i2ndstr.gt.0) call secstrp2dihc
1326 c      call geom_to_var(nvar,x)
1327 c      call etotal(energia(0))
1328 c      call enerprint(energia(0))
1329 c      call briefout(0,etot)
1330 c      stop
1331 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1332 cd    write (iout,'(a)') 'Variable list:'
1333 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1334 #ifdef MPI
1335       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1336      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1337      &  'Processor',myrank,': end reading molecular data.'
1338 #endif
1339       return
1340       end
1341 c--------------------------------------------------------------------------
1342       logical function seq_comp(itypea,itypeb,length)
1343       implicit none
1344       integer length,itypea(length),itypeb(length)
1345       integer i
1346       do i=1,length
1347         if (itypea(i).ne.itypeb(i)) then
1348           seq_comp=.false.
1349           return
1350         endif
1351       enddo
1352       seq_comp=.true.
1353       return
1354       end
1355 c-----------------------------------------------------------------------------
1356       subroutine read_bridge
1357 C Read information about disulfide bridges.
1358       implicit real*8 (a-h,o-z)
1359       include 'DIMENSIONS'
1360 #ifdef MPI
1361       include 'mpif.h'
1362 #endif
1363       include 'COMMON.IOUNITS'
1364       include 'COMMON.GEO'
1365       include 'COMMON.VAR'
1366       include 'COMMON.INTERACT'
1367       include 'COMMON.LOCAL'
1368       include 'COMMON.NAMES'
1369       include 'COMMON.CHAIN'
1370       include 'COMMON.FFIELD'
1371       include 'COMMON.SBRIDGE'
1372       include 'COMMON.HEADER'
1373       include 'COMMON.CONTROL'
1374       include 'COMMON.DBASE'
1375       include 'COMMON.THREAD'
1376       include 'COMMON.TIME1'
1377       include 'COMMON.SETUP'
1378 C Read bridging residues.
1379       read (inp,*) ns,(iss(i),i=1,ns)
1380       print *,'ns=',ns
1381       if(me.eq.king.or..not.out1file)
1382      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1383 C Check whether the specified bridging residues are cystines.
1384       do i=1,ns
1385         if (itype(iss(i)).ne.1) then
1386           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1387      &   'Do you REALLY think that the residue ',
1388      &    restyp(itype(iss(i))),i,
1389      &   ' can form a disulfide bridge?!!!'
1390           write (*,'(2a,i3,a)')
1391      &   'Do you REALLY think that the residue ',
1392      &    restyp(itype(iss(i))),i,
1393      &   ' can form a disulfide bridge?!!!'
1394 #ifdef MPI
1395          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1396          stop
1397 #endif
1398         endif
1399       enddo
1400 C Read preformed bridges.
1401       if (ns.gt.0) then
1402       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1403       if(fg_rank.eq.0)
1404      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1405       if (nss.gt.0) then
1406         nhpb=nss
1407 C Check if the residues involved in bridges are in the specified list of
1408 C bridging residues.
1409         do i=1,nss
1410           do j=1,i-1
1411             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1412      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1413               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1414      &      ' contains residues present in other pairs.'
1415               write (*,'(a,i3,a)') 'Disulfide pair',i,
1416      &      ' contains residues present in other pairs.'
1417 #ifdef MPI
1418               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1419               stop
1420 #endif
1421             endif
1422           enddo
1423           do j=1,ns
1424             if (ihpb(i).eq.iss(j)) goto 10
1425           enddo
1426           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1427    10     continue
1428           do j=1,ns
1429             if (jhpb(i).eq.iss(j)) goto 20
1430           enddo
1431           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1432    20     continue
1433           dhpb(i)=dbr
1434           forcon(i)=fbr
1435         enddo
1436         do i=1,nss
1437           ihpb(i)=ihpb(i)+nres
1438           jhpb(i)=jhpb(i)+nres
1439         enddo
1440       endif
1441       endif
1442       return
1443       end
1444 c----------------------------------------------------------------------------
1445       subroutine read_x(kanal,*)
1446       implicit real*8 (a-h,o-z)
1447       include 'DIMENSIONS'
1448       include 'COMMON.GEO'
1449       include 'COMMON.VAR'
1450       include 'COMMON.CHAIN'
1451       include 'COMMON.IOUNITS'
1452       include 'COMMON.CONTROL'
1453       include 'COMMON.LOCAL'
1454       include 'COMMON.INTERACT'
1455 c Read coordinates from input
1456 c
1457       read(kanal,'(8f10.5)',end=10,err=10)
1458      &  ((c(l,k),l=1,3),k=1,nres),
1459      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1460       do j=1,3
1461         c(j,nres+1)=c(j,1)
1462         c(j,2*nres)=c(j,nres)
1463       enddo
1464       call int_from_cart1(.false.)
1465       do i=1,nres-1
1466         do j=1,3
1467           dc(j,i)=c(j,i+1)-c(j,i)
1468           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1469         enddo
1470       enddo
1471       do i=nnt,nct
1472         if (itype(i).ne.10) then
1473           do j=1,3
1474             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1475             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1476           enddo
1477         endif
1478       enddo
1479
1480       return
1481    10 return1
1482       end
1483 c----------------------------------------------------------------------------
1484       subroutine read_threadbase
1485       implicit real*8 (a-h,o-z)
1486       include 'DIMENSIONS'
1487       include 'COMMON.IOUNITS'
1488       include 'COMMON.GEO'
1489       include 'COMMON.VAR'
1490       include 'COMMON.INTERACT'
1491       include 'COMMON.LOCAL'
1492       include 'COMMON.NAMES'
1493       include 'COMMON.CHAIN'
1494       include 'COMMON.FFIELD'
1495       include 'COMMON.SBRIDGE'
1496       include 'COMMON.HEADER'
1497       include 'COMMON.CONTROL'
1498       include 'COMMON.DBASE'
1499       include 'COMMON.THREAD'
1500       include 'COMMON.TIME1'
1501 C Read pattern database for threading.
1502       read (icbase,*) nseq
1503       do i=1,nseq
1504         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1505      &   nres_base(2,i),nres_base(3,i)
1506         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1507      &   nres_base(1,i))
1508 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1509 c    &   nres_base(2,i),nres_base(3,i)
1510 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1511 c    &   nres_base(1,i))
1512       enddo
1513       close (icbase)
1514       if (weidis.eq.0.0D0) weidis=0.1D0
1515       do i=nnt,nct
1516         do j=i+2,nct
1517           nhpb=nhpb+1
1518           ihpb(nhpb)=i
1519           jhpb(nhpb)=j
1520           forcon(nhpb)=weidis
1521         enddo
1522       enddo
1523       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1524       write (iout,'(a,i5)') 'nexcl: ',nexcl
1525       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1526       return
1527       end
1528 c------------------------------------------------------------------------------
1529       subroutine setup_var
1530       implicit real*8 (a-h,o-z)
1531       include 'DIMENSIONS'
1532       include 'COMMON.IOUNITS'
1533       include 'COMMON.GEO'
1534       include 'COMMON.VAR'
1535       include 'COMMON.INTERACT'
1536       include 'COMMON.LOCAL'
1537       include 'COMMON.NAMES'
1538       include 'COMMON.CHAIN'
1539       include 'COMMON.FFIELD'
1540       include 'COMMON.SBRIDGE'
1541       include 'COMMON.HEADER'
1542       include 'COMMON.CONTROL'
1543       include 'COMMON.DBASE'
1544       include 'COMMON.THREAD'
1545       include 'COMMON.TIME1'
1546 C Set up variable list.
1547       ntheta=nres-2
1548       nphi=nres-3
1549       nvar=ntheta+nphi
1550       nside=0
1551       do i=2,nres-1
1552         if (itype(i).ne.10) then
1553           nside=nside+1
1554           ialph(i,1)=nvar+nside
1555           ialph(nside,2)=i
1556         endif
1557       enddo
1558       if (indphi.gt.0) then
1559         nvar=nphi
1560       else if (indback.gt.0) then
1561         nvar=nphi+ntheta
1562       else
1563         nvar=nvar+2*nside
1564       endif
1565 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1566       return
1567       end
1568 c----------------------------------------------------------------------------
1569       subroutine gen_dist_constr
1570 C Generate CA distance constraints.
1571       implicit real*8 (a-h,o-z)
1572       include 'DIMENSIONS'
1573       include 'COMMON.IOUNITS'
1574       include 'COMMON.GEO'
1575       include 'COMMON.VAR'
1576       include 'COMMON.INTERACT'
1577       include 'COMMON.LOCAL'
1578       include 'COMMON.NAMES'
1579       include 'COMMON.CHAIN'
1580       include 'COMMON.FFIELD'
1581       include 'COMMON.SBRIDGE'
1582       include 'COMMON.HEADER'
1583       include 'COMMON.CONTROL'
1584       include 'COMMON.DBASE'
1585       include 'COMMON.THREAD'
1586       include 'COMMON.TIME1'
1587       dimension itype_pdb(maxres)
1588       common /pizda/ itype_pdb
1589       character*2 iden
1590 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1591 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1592 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1593 cd     & ' nsup',nsup
1594       do i=nstart_sup,nstart_sup+nsup-1
1595 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1596 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1597         do j=i+2,nstart_sup+nsup-1
1598           nhpb=nhpb+1
1599           ihpb(nhpb)=i+nstart_seq-nstart_sup
1600           jhpb(nhpb)=j+nstart_seq-nstart_sup
1601           forcon(nhpb)=weidis
1602           dhpb(nhpb)=dist(i,j)
1603         enddo
1604       enddo
1605 cd      write (iout,'(a)') 'Distance constraints:'
1606 cd      do i=nss+1,nhpb
1607 cd        ii=ihpb(i)
1608 cd        jj=jhpb(i)
1609 cd        iden='CA'
1610 cd        if (ii.gt.nres) then
1611 cd          iden='SC'
1612 cd          ii=ii-nres
1613 cd          jj=jj-nres
1614 cd        endif
1615 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1616 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1617 cd     &  dhpb(i),forcon(i)
1618 cd      enddo
1619       return
1620       end
1621 c----------------------------------------------------------------------------
1622       subroutine map_read
1623       implicit real*8 (a-h,o-z)
1624       include 'DIMENSIONS'
1625       include 'COMMON.MAP'
1626       include 'COMMON.IOUNITS'
1627       character*3 angid(4) /'THE','PHI','ALP','OME'/
1628       character*80 mapcard,ucase
1629       do imap=1,nmap
1630         read (inp,'(a)') mapcard
1631         mapcard=ucase(mapcard)
1632         if (index(mapcard,'PHI').gt.0) then
1633           kang(imap)=1
1634         else if (index(mapcard,'THE').gt.0) then
1635           kang(imap)=2
1636         else if (index(mapcard,'ALP').gt.0) then
1637           kang(imap)=3
1638         else if (index(mapcard,'OME').gt.0) then
1639           kang(imap)=4
1640         else
1641           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1642           stop 'Error - illegal variable spec in MAP card.'
1643         endif
1644         call readi (mapcard,'RES1',res1(imap),0)
1645         call readi (mapcard,'RES2',res2(imap),0)
1646         if (res1(imap).eq.0) then
1647           res1(imap)=res2(imap)
1648         else if (res2(imap).eq.0) then
1649           res2(imap)=res1(imap)
1650         endif
1651         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1652           write (iout,'(a)')
1653      &    'Error - illegal definition of variable group in MAP.'
1654           stop 'Error - illegal definition of variable group in MAP.'
1655         endif
1656         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1657         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1658         call readi(mapcard,'NSTEP',nstep(imap),0)
1659         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1660           write (iout,'(a)')
1661      &     'Illegal boundary and/or step size specification in MAP.'
1662           stop 'Illegal boundary and/or step size specification in MAP.'
1663         endif
1664       enddo ! imap
1665       return
1666       end
1667 c----------------------------------------------------------------------------
1668 csa      subroutine csaread
1669 csa      implicit real*8 (a-h,o-z)
1670 csa      include 'DIMENSIONS'
1671 csa      include 'COMMON.IOUNITS'
1672 csa      include 'COMMON.GEO'
1673 csa      include 'COMMON.CSA'
1674 csa      include 'COMMON.BANK'
1675 csa      include 'COMMON.CONTROL'
1676 csa      character*80 ucase
1677 csa      character*620 mcmcard
1678 csa      call card_concat(mcmcard)
1679 csa
1680 csa      call readi(mcmcard,'NCONF',nconf,50)
1681 csa      call readi(mcmcard,'NADD',nadd,0)
1682 csa      call readi(mcmcard,'JSTART',jstart,1)
1683 csa      call readi(mcmcard,'JEND',jend,1)
1684 csa      call readi(mcmcard,'NSTMAX',nstmax,500000)
1685 csa      call readi(mcmcard,'N0',n0,1)
1686 csa      call readi(mcmcard,'N1',n1,6)
1687 csa      call readi(mcmcard,'N2',n2,4)
1688 csa      call readi(mcmcard,'N3',n3,0)
1689 csa      call readi(mcmcard,'N4',n4,0)
1690 csa      call readi(mcmcard,'N5',n5,0)
1691 csa      call readi(mcmcard,'N6',n6,10)
1692 csa      call readi(mcmcard,'N7',n7,0)
1693 csa      call readi(mcmcard,'N8',n8,0)
1694 csa      call readi(mcmcard,'N9',n9,0)
1695 csa      call readi(mcmcard,'N14',n14,0)
1696 csa      call readi(mcmcard,'N15',n15,0)
1697 csa      call readi(mcmcard,'N16',n16,0)
1698 csa      call readi(mcmcard,'N17',n17,0)
1699 csa      call readi(mcmcard,'N18',n18,0)
1700 csa
1701 csa      vdisulf=(index(mcmcard,'DYNSS').gt.0)
1702 csa
1703 csa      call readi(mcmcard,'NDIFF',ndiff,2)
1704 csa      call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1705 csa      call readi(mcmcard,'IS1',is1,1)
1706 csa      call readi(mcmcard,'IS2',is2,8)
1707 csa      call readi(mcmcard,'NRAN0',nran0,4)
1708 csa      call readi(mcmcard,'NRAN1',nran1,2)
1709 csa      call readi(mcmcard,'IRR',irr,1)
1710 csa      call readi(mcmcard,'NSEED',nseed,20)
1711 csa      call readi(mcmcard,'NTOTAL',ntotal,10000)
1712 csa      call reada(mcmcard,'CUT1',cut1,2.0d0)
1713 csa      call reada(mcmcard,'CUT2',cut2,5.0d0)
1714 csa      call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1715 csa      call readi(mcmcard,'ICMAX',icmax,3)
1716 csa      call readi(mcmcard,'IRESTART',irestart,0)
1717 csac!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1718 csa      ntbankm=0
1719 csac!bankt
1720 csa      call reada(mcmcard,'DELE',dele,20.0d0)
1721 csa      call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1722 csa      call readi(mcmcard,'IREF',iref,0)
1723 csa      call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1724 csa      call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1725 csa      call readi(mcmcard,'NCONF_IN',nconf_in,0)
1726 csa      call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1727 csa      write (iout,*) "NCONF_IN",nconf_in
1728 csa      return
1729 csa      end
1730 c----------------------------------------------------------------------------
1731 cfmc      subroutine mcmfread
1732 cfmc      implicit real*8 (a-h,o-z)
1733 cfmc      include 'DIMENSIONS'
1734 cfmc      include 'COMMON.MCMF'
1735 cfmc      include 'COMMON.IOUNITS'
1736 cfmc      include 'COMMON.GEO'
1737 cfmc      character*80 ucase
1738 cfmc      character*620 mcmcard
1739 cfmc      call card_concat(mcmcard)
1740 cfmc
1741 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1742 cfmc      write(iout,*)'MAXRANT=',maxrant
1743 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1744 cfmc      write(iout,*)'MAXFAM=',maxfam
1745 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1746 cfmc      write(iout,*)'NNET1=',nnet1
1747 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1748 cfmc      write(iout,*)'NNET2=',nnet2
1749 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1750 cfmc      write(iout,*)'NNET3=',nnet3
1751 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1752 cfmc      write(iout,*)'ILASTT=',ilastt
1753 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1754 cfmc      write(iout,*)'MAXSTR=',maxstr
1755 cfmc      maxstr_f=maxstr/maxfam
1756 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1757 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1758 cfmc      write(iout,*)'NMCMF=',nmcmf
1759 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1760 cfmc      write(iout,*)'IFOCUS=',ifocus
1761 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1762 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1763 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1764 cfmc      write(iout,*)'INTPRT=',intprt
1765 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1766 cfmc      write(iout,*)'IPRT=',iprt
1767 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1768 cfmc      write(iout,*)'IMAXTR=',imaxtr
1769 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1770 cfmc      write(iout,*)'MAXEVEN=',maxeven
1771 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1772 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1773 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1774 cfmc      write(iout,*)'INIMIN=',inimin
1775 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1776 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1777 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1778 cfmc      write(iout,*)'NTHREAD=',nthread
1779 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1780 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1781 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1782 cfmc      write(iout,*)'MAXPERT=',maxpert
1783 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1784 cfmc      write(iout,*)'IRMSD=',irmsd
1785 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1786 cfmc      write(iout,*)'DENEMIN=',denemin
1787 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1788 cfmc      write(iout,*)'RCUT1S=',rcut1s
1789 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1790 cfmc      write(iout,*)'RCUT1E=',rcut1e
1791 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1792 cfmc      write(iout,*)'RCUT2S=',rcut2s
1793 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1794 cfmc      write(iout,*)'RCUT2E=',rcut2e
1795 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1796 cfmc      write(iout,*)'DPERT1=',d_pert1
1797 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1798 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1799 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1800 cfmc      write(iout,*)'DPERT2=',d_pert2
1801 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1802 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1803 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1804 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1805 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1806 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1807 cfmc      d_pert1=deg2rad*d_pert1
1808 cfmc      d_pert1a=deg2rad*d_pert1a
1809 cfmc      d_pert2=deg2rad*d_pert2
1810 cfmc      d_pert2a=deg2rad*d_pert2a
1811 cfmc      d_pert2b=deg2rad*d_pert2b
1812 cfmc      d_pert2c=deg2rad*d_pert2c
1813 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1814 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1815 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1816 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1817 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1818 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1819 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1820 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1821 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1822 cfmc      write(iout,*)'RCUTINI=',rcutini
1823 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1824 cfmc      write(iout,*)'GRAT=',grat
1825 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1826 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1827 cfmc
1828 cfmc      return
1829 cfmc      end
1830 c----------------------------------------------------------------------------
1831       subroutine mcmread
1832       implicit real*8 (a-h,o-z)
1833       include 'DIMENSIONS'
1834       include 'COMMON.MCM'
1835       include 'COMMON.MCE'
1836       include 'COMMON.IOUNITS'
1837       character*80 ucase
1838       character*320 mcmcard
1839       call card_concat(mcmcard)
1840       call readi(mcmcard,'MAXACC',maxacc,100)
1841       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1842       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1843       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1844       call readi(mcmcard,'MAXREPM',maxrepm,200)
1845       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1846       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1847       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1848       call reada(mcmcard,'E_UP',e_up,5.0D0)
1849       call reada(mcmcard,'DELTE',delte,0.1D0)
1850       call readi(mcmcard,'NSWEEP',nsweep,5)
1851       call readi(mcmcard,'NSTEPH',nsteph,0)
1852       call readi(mcmcard,'NSTEPC',nstepc,0)
1853       call reada(mcmcard,'TMIN',tmin,298.0D0)
1854       call reada(mcmcard,'TMAX',tmax,298.0D0)
1855       call readi(mcmcard,'NWINDOW',nwindow,0)
1856       call readi(mcmcard,'PRINT_MC',print_mc,0)
1857       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1858       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1859       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1860       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1861       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1862       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1863       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1864       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1865       if (nwindow.gt.0) then
1866         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1867         do i=1,nwindow
1868           winlen(i)=winend(i)-winstart(i)+1
1869         enddo
1870       endif
1871       if (tmax.lt.tmin) tmax=tmin
1872       if (tmax.eq.tmin) then
1873         nstepc=0
1874         nsteph=0
1875       endif
1876       if (nstepc.gt.0 .and. nsteph.gt.0) then
1877         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1878         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1879       endif
1880 C Probabilities of different move types
1881       sumpro_type(0)=0.0D0
1882       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1883       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1884       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1885       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1886       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1887       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1888       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1889       do i=1,MaxMoveType
1890         print *,'i',i,' sumprotype',sumpro_type(i)
1891         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1892         print *,'i',i,' sumprotype',sumpro_type(i)
1893       enddo
1894       return
1895       end
1896 c----------------------------------------------------------------------------
1897       subroutine read_minim
1898       implicit real*8 (a-h,o-z)
1899       include 'DIMENSIONS'
1900       include 'COMMON.MINIM'
1901       include 'COMMON.IOUNITS'
1902       character*80 ucase
1903       character*320 minimcard
1904       call card_concat(minimcard)
1905       call readi(minimcard,'MAXMIN',maxmin,2000)
1906       call readi(minimcard,'MAXFUN',maxfun,5000)
1907       call readi(minimcard,'MINMIN',minmin,maxmin)
1908       call readi(minimcard,'MINFUN',minfun,maxmin)
1909       call reada(minimcard,'TOLF',tolf,1.0D-2)
1910       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1911       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1912       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1913       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1914       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1915      &         'Options in energy minimization:'
1916       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1917      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1918      & 'MinMin:',MinMin,' MinFun:',MinFun,
1919      & ' TolF:',TolF,' RTolF:',RTolF
1920       return
1921       end
1922 c----------------------------------------------------------------------------
1923       subroutine read_angles(kanal,*)
1924       implicit real*8 (a-h,o-z)
1925       include 'DIMENSIONS'
1926       include 'COMMON.GEO'
1927       include 'COMMON.VAR'
1928       include 'COMMON.CHAIN'
1929       include 'COMMON.IOUNITS'
1930       include 'COMMON.CONTROL'
1931 c Read angles from input
1932 c
1933        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1934        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1935        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1936        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1937
1938        do i=1,nres
1939 c 9/7/01 avoid 180 deg valence angle
1940         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1941 c
1942         theta(i)=deg2rad*theta(i)
1943         phi(i)=deg2rad*phi(i)
1944         alph(i)=deg2rad*alph(i)
1945         omeg(i)=deg2rad*omeg(i)
1946        enddo
1947       return
1948    10 return1
1949       end
1950 c----------------------------------------------------------------------------
1951       subroutine reada(rekord,lancuch,wartosc,default)
1952       implicit none
1953       character*(*) rekord,lancuch
1954       double precision wartosc,default
1955       integer ilen,iread
1956       external ilen
1957       iread=index(rekord,lancuch)
1958       if (iread.eq.0) then
1959         wartosc=default
1960         return
1961       endif
1962       iread=iread+ilen(lancuch)+1
1963       read (rekord(iread:),*,err=10,end=10) wartosc
1964       return
1965   10  wartosc=default
1966       return
1967       end
1968 c----------------------------------------------------------------------------
1969       subroutine readi(rekord,lancuch,wartosc,default)
1970       implicit none
1971       character*(*) rekord,lancuch
1972       integer wartosc,default
1973       integer ilen,iread
1974       external ilen
1975       iread=index(rekord,lancuch)
1976       if (iread.eq.0) then
1977         wartosc=default
1978         return
1979       endif
1980       iread=iread+ilen(lancuch)+1
1981       read (rekord(iread:),*,err=10,end=10) wartosc
1982       return
1983   10  wartosc=default
1984       return
1985       end
1986 c----------------------------------------------------------------------------
1987       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1988       implicit none
1989       integer dim,i
1990       integer tablica(dim),default
1991       character*(*) rekord,lancuch
1992       character*80 aux
1993       integer ilen,iread
1994       external ilen
1995       do i=1,dim
1996         tablica(i)=default
1997       enddo
1998       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1999       if (iread.eq.0) return
2000       iread=iread+ilen(lancuch)+1
2001       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2002    10 return
2003       end
2004 c----------------------------------------------------------------------------
2005       subroutine multreada(rekord,lancuch,tablica,dim,default)
2006       implicit none
2007       integer dim,i
2008       double precision tablica(dim),default
2009       character*(*) rekord,lancuch
2010       character*80 aux
2011       integer ilen,iread
2012       external ilen
2013       do i=1,dim
2014         tablica(i)=default
2015       enddo
2016       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2017       if (iread.eq.0) return
2018       iread=iread+ilen(lancuch)+1
2019       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2020    10 return
2021       end
2022 c----------------------------------------------------------------------------
2023       subroutine openunits
2024       implicit real*8 (a-h,o-z)
2025       include 'DIMENSIONS'
2026 #ifdef MPI
2027       include 'mpif.h'
2028       character*16 form,nodename
2029       integer nodelen
2030 #endif
2031       include 'COMMON.SETUP'
2032       include 'COMMON.IOUNITS'
2033       include 'COMMON.MD'
2034       include 'COMMON.CONTROL'
2035       integer lenpre,lenpot,ilen,lentmp
2036       external ilen
2037       character*3 out1file_text,ucase
2038       character*3 ll
2039       external ucase
2040 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2041       call getenv_loc("PREFIX",prefix)
2042       pref_orig = prefix
2043       call getenv_loc("POT",pot)
2044       call getenv_loc("DIRTMP",tmpdir)
2045       call getenv_loc("CURDIR",curdir)
2046       call getenv_loc("OUT1FILE",out1file_text)
2047 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2048       out1file_text=ucase(out1file_text)
2049       if (out1file_text(1:1).eq."Y") then
2050         out1file=.true.
2051       else
2052         out1file=fg_rank.gt.0
2053       endif
2054       lenpre=ilen(prefix)
2055       lenpot=ilen(pot)
2056       lentmp=ilen(tmpdir)
2057       if (lentmp.gt.0) then
2058           write (*,'(80(1h!))')
2059           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2060           write (*,'(80(1h!))')
2061           write (*,*)"All output files will be on node /tmp directory."
2062 #ifdef MPI
2063         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2064         if (me.eq.king) then
2065           write (*,*) "The master node is ",nodename
2066         else if (fg_rank.eq.0) then
2067           write (*,*) "I am the CG slave node ",nodename
2068         else
2069           write (*,*) "I am the FG slave node ",nodename
2070         endif
2071 #endif
2072         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2073         lenpre = lentmp+lenpre+1
2074       endif
2075       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2076 C Get the names and open the input files
2077 #if defined(WINIFL) || defined(WINPGI)
2078       open(1,file=pref_orig(:ilen(pref_orig))//
2079      &  '.inp',status='old',readonly,shared)
2080        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2081 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2082 C Get parameter filenames and open the parameter files.
2083       call getenv_loc('BONDPAR',bondname)
2084       open (ibond,file=bondname,status='old',readonly,shared)
2085       call getenv_loc('THETPAR',thetname)
2086       open (ithep,file=thetname,status='old',readonly,shared)
2087 #ifndef CRYST_THETA
2088       call getenv_loc('THETPARPDB',thetname_pdb)
2089       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2090 #endif
2091       call getenv_loc('ROTPAR',rotname)
2092       open (irotam,file=rotname,status='old',readonly,shared)
2093 #ifndef CRYST_SC
2094       call getenv_loc('ROTPARPDB',rotname_pdb)
2095       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2096 #endif
2097       call getenv_loc('TORPAR',torname)
2098       open (itorp,file=torname,status='old',readonly,shared)
2099       call getenv_loc('TORDPAR',tordname)
2100       open (itordp,file=tordname,status='old',readonly,shared)
2101       call getenv_loc('FOURIER',fouriername)
2102       open (ifourier,file=fouriername,status='old',readonly,shared)
2103       call getenv_loc('ELEPAR',elename)
2104       open (ielep,file=elename,status='old',readonly,shared)
2105       call getenv_loc('SIDEPAR',sidename)
2106       open (isidep,file=sidename,status='old',readonly,shared)
2107 #elif (defined CRAY) || (defined AIX)
2108       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2109      &  action='read')
2110 c      print *,"Processor",myrank," opened file 1"
2111       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2112 c      print *,"Processor",myrank," opened file 9"
2113 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2114 C Get parameter filenames and open the parameter files.
2115       call getenv_loc('BONDPAR',bondname)
2116       open (ibond,file=bondname,status='old',action='read')
2117 c      print *,"Processor",myrank," opened file IBOND"
2118       call getenv_loc('THETPAR',thetname)
2119       open (ithep,file=thetname,status='old',action='read')
2120 c      print *,"Processor",myrank," opened file ITHEP"
2121 #ifndef CRYST_THETA
2122       call getenv_loc('THETPARPDB',thetname_pdb)
2123       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2124 #endif
2125       call getenv_loc('ROTPAR',rotname)
2126       open (irotam,file=rotname,status='old',action='read')
2127 c      print *,"Processor",myrank," opened file IROTAM"
2128 #ifndef CRYST_SC
2129       call getenv_loc('ROTPARPDB',rotname_pdb)
2130       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2131 #endif
2132       call getenv_loc('TORPAR',torname)
2133       open (itorp,file=torname,status='old',action='read')
2134 c      print *,"Processor",myrank," opened file ITORP"
2135       call getenv_loc('TORDPAR',tordname)
2136       open (itordp,file=tordname,status='old',action='read')
2137 c      print *,"Processor",myrank," opened file ITORDP"
2138       call getenv_loc('SCCORPAR',sccorname)
2139       open (isccor,file=sccorname,status='old',action='read')
2140 c      print *,"Processor",myrank," opened file ISCCOR"
2141       call getenv_loc('FOURIER',fouriername)
2142       open (ifourier,file=fouriername,status='old',action='read')
2143 c      print *,"Processor",myrank," opened file IFOURIER"
2144       call getenv_loc('ELEPAR',elename)
2145       open (ielep,file=elename,status='old',action='read')
2146 c      print *,"Processor",myrank," opened file IELEP"
2147       call getenv_loc('SIDEPAR',sidename)
2148       open (isidep,file=sidename,status='old',action='read')
2149 c      print *,"Processor",myrank," opened file ISIDEP"
2150 c      print *,"Processor",myrank," opened parameter files"
2151 #elif (defined G77)
2152       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2153       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2154 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2155 C Get parameter filenames and open the parameter files.
2156       call getenv_loc('BONDPAR',bondname)
2157       open (ibond,file=bondname,status='old')
2158       call getenv_loc('THETPAR',thetname)
2159       open (ithep,file=thetname,status='old')
2160 #ifndef CRYST_THETA
2161       call getenv_loc('THETPARPDB',thetname_pdb)
2162       open (ithep_pdb,file=thetname_pdb,status='old')
2163 #endif
2164       call getenv_loc('ROTPAR',rotname)
2165       open (irotam,file=rotname,status='old')
2166 #ifndef CRYST_SC
2167       call getenv_loc('ROTPARPDB',rotname_pdb)
2168       open (irotam_pdb,file=rotname_pdb,status='old')
2169 #endif
2170       call getenv_loc('TORPAR',torname)
2171       open (itorp,file=torname,status='old')
2172       call getenv_loc('TORDPAR',tordname)
2173       open (itordp,file=tordname,status='old')
2174       call getenv_loc('SCCORPAR',sccorname)
2175       open (isccor,file=sccorname,status='old')
2176       call getenv_loc('FOURIER',fouriername)
2177       open (ifourier,file=fouriername,status='old')
2178       call getenv_loc('ELEPAR',elename)
2179       open (ielep,file=elename,status='old')
2180       call getenv_loc('SIDEPAR',sidename)
2181       open (isidep,file=sidename,status='old')
2182 #else
2183       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2184      &action='read')
2185        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2186 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2187 C Get parameter filenames and open the parameter files.
2188       call getenv_loc('BONDPAR',bondname)
2189       open (ibond,file=bondname,status='old',action='read')
2190       call getenv_loc('THETPAR',thetname)
2191       open (ithep,file=thetname,status='old',action='read')
2192 #ifndef CRYST_THETA
2193       call getenv_loc('THETPARPDB',thetname_pdb)
2194       print *,"thetname_pdb ",thetname_pdb
2195       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2196       print *,ithep_pdb," opened"
2197 #endif
2198       call getenv_loc('ROTPAR',rotname)
2199       open (irotam,file=rotname,status='old',action='read')
2200 #ifndef CRYST_SC
2201       call getenv_loc('ROTPARPDB',rotname_pdb)
2202       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2203 #endif
2204       call getenv_loc('TORPAR',torname)
2205       open (itorp,file=torname,status='old',action='read')
2206       call getenv_loc('TORDPAR',tordname)
2207       open (itordp,file=tordname,status='old',action='read')
2208       call getenv_loc('SCCORPAR',sccorname)
2209       open (isccor,file=sccorname,status='old',action='read')
2210       call getenv_loc('FOURIER',fouriername)
2211       open (ifourier,file=fouriername,status='old',action='read')
2212       call getenv_loc('ELEPAR',elename)
2213       open (ielep,file=elename,status='old',action='read')
2214       call getenv_loc('SIDEPAR',sidename)
2215       open (isidep,file=sidename,status='old',action='read')
2216 #endif
2217 #ifndef OLDSCP
2218 C
2219 C 8/9/01 In the newest version SCp interaction constants are read from a file
2220 C Use -DOLDSCP to use hard-coded constants instead.
2221 C
2222       call getenv_loc('SCPPAR',scpname)
2223 #if defined(WINIFL) || defined(WINPGI)
2224       open (iscpp,file=scpname,status='old',readonly,shared)
2225 #elif (defined CRAY)  || (defined AIX)
2226       open (iscpp,file=scpname,status='old',action='read')
2227 #elif (defined G77)
2228       open (iscpp,file=scpname,status='old')
2229 #else
2230       open (iscpp,file=scpname,status='old',action='read')
2231 #endif
2232 #endif
2233       call getenv_loc('PATTERN',patname)
2234 #if defined(WINIFL) || defined(WINPGI)
2235       open (icbase,file=patname,status='old',readonly,shared)
2236 #elif (defined CRAY)  || (defined AIX)
2237       open (icbase,file=patname,status='old',action='read')
2238 #elif (defined G77)
2239       open (icbase,file=patname,status='old')
2240 #else
2241       open (icbase,file=patname,status='old',action='read')
2242 #endif
2243 #ifdef MPI
2244 C Open output file only for CG processes
2245 c      print *,"Processor",myrank," fg_rank",fg_rank
2246       if (fg_rank.eq.0) then
2247
2248       if (nodes.eq.1) then
2249         npos=3
2250       else
2251         npos = dlog10(dfloat(nodes-1))+1
2252       endif
2253       if (npos.lt.3) npos=3
2254       write (liczba,'(i1)') npos
2255       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2256      &  //')'
2257       write (liczba,form) me
2258       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2259      &  liczba(:ilen(liczba))
2260       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2261      &  //'.int'
2262       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2263      &  //'.pdb'
2264       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2265      &  liczba(:ilen(liczba))//'.mol2'
2266       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2267      &  liczba(:ilen(liczba))//'.stat'
2268       if (lentmp.gt.0)
2269      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2270      &      //liczba(:ilen(liczba))//'.stat')
2271       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2272      &  //'.rst'
2273       if(usampl) then
2274           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2275      & liczba(:ilen(liczba))//'.const'
2276       endif
2277
2278       endif
2279 #else
2280       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2281       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2282       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2283       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2284       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2285       if (lentmp.gt.0)
2286      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2287      &    //'.stat')
2288       rest2name=prefix(:ilen(prefix))//'.rst'
2289       if(usampl) then
2290          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2291       endif
2292 #endif
2293 #if defined(AIX) || defined(PGI)
2294       if (me.eq.king .or. .not. out1file)
2295      &   open(iout,file=outname,status='unknown')
2296 c#define DEBUG
2297 #ifdef DEBUG
2298       if (fg_rank.gt.0) then
2299         write (liczba,'(i3.3)') myrank/nfgtasks
2300         write (ll,'(bz,i3.3)') fg_rank
2301         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2302      &   status='unknown')
2303       endif
2304 #endif
2305 c#undef DEBUG
2306       if(me.eq.king) then
2307        open(igeom,file=intname,status='unknown',position='append')
2308        open(ipdb,file=pdbname,status='unknown')
2309        open(imol2,file=mol2name,status='unknown')
2310        open(istat,file=statname,status='unknown',position='append')
2311       else
2312 c1out       open(iout,file=outname,status='unknown')
2313       endif
2314 #else
2315       if (me.eq.king .or. .not.out1file)
2316      &    open(iout,file=outname,status='unknown')
2317 c#define DEBUG
2318 #ifdef DEBUG
2319       if (fg_rank.gt.0) then
2320         print "Processor",fg_rank," opening output file"
2321         write (liczba,'(i3.3)') myrank/nfgtasks
2322         write (ll,'(bz,i3.3)') fg_rank
2323         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2324      &   status='unknown')
2325       endif
2326 #endif
2327 c#undef DEBUG
2328       if(me.eq.king) then
2329        open(igeom,file=intname,status='unknown',access='append')
2330        open(ipdb,file=pdbname,status='unknown')
2331        open(imol2,file=mol2name,status='unknown')
2332        open(istat,file=statname,status='unknown',access='append')
2333       else
2334 c1out       open(iout,file=outname,status='unknown')
2335       endif
2336 #endif
2337 csa      csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2338 csa      csa_seed=prefix(:lenpre)//'.CSA.seed'
2339 csa      csa_history=prefix(:lenpre)//'.CSA.history'
2340 csa      csa_bank=prefix(:lenpre)//'.CSA.bank'
2341 csa      csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2342 csa      csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2343 csa      csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2344 csac!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2345 csa      csa_int=prefix(:lenpre)//'.int'
2346 csa      csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2347 csa      csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2348 csa      csa_in=prefix(:lenpre)//'.CSA.in'
2349 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2350 C Write file names
2351       if (me.eq.king)then
2352       write (iout,'(80(1h-))')
2353       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2354       write (iout,'(80(1h-))')
2355       write (iout,*) "Input file                      : ",
2356      &  pref_orig(:ilen(pref_orig))//'.inp'
2357       write (iout,*) "Output file                     : ",
2358      &  outname(:ilen(outname))
2359       write (iout,*)
2360       write (iout,*) "Sidechain potential file        : ",
2361      &  sidename(:ilen(sidename))
2362 #ifndef OLDSCP
2363       write (iout,*) "SCp potential file              : ",
2364      &  scpname(:ilen(scpname))
2365 #endif
2366       write (iout,*) "Electrostatic potential file    : ",
2367      &  elename(:ilen(elename))
2368       write (iout,*) "Cumulant coefficient file       : ",
2369      &  fouriername(:ilen(fouriername))
2370       write (iout,*) "Torsional parameter file        : ",
2371      &  torname(:ilen(torname))
2372       write (iout,*) "Double torsional parameter file : ",
2373      &  tordname(:ilen(tordname))
2374       write (iout,*) "SCCOR parameter file : ",
2375      &  sccorname(:ilen(sccorname))
2376       write (iout,*) "Bond & inertia constant file    : ",
2377      &  bondname(:ilen(bondname))
2378       write (iout,*) "Bending parameter file          : ",
2379      &  thetname(:ilen(thetname))
2380       write (iout,*) "Rotamer parameter file          : ",
2381      &  rotname(:ilen(rotname))
2382       write (iout,*) "Threading database              : ",
2383      &  patname(:ilen(patname))
2384       if (lentmp.ne.0)
2385      &write (iout,*)" DIRTMP                          : ",
2386      &  tmpdir(:lentmp)
2387       write (iout,'(80(1h-))')
2388       endif
2389       return
2390       end
2391 c----------------------------------------------------------------------------
2392       subroutine card_concat(card)
2393       implicit real*8 (a-h,o-z)
2394       include 'DIMENSIONS'
2395       include 'COMMON.IOUNITS'
2396       character*(*) card
2397       character*80 karta,ucase
2398       external ilen
2399       read (inp,'(a)') karta
2400       karta=ucase(karta)
2401       card=' '
2402       do while (karta(80:80).eq.'&')
2403         card=card(:ilen(card)+1)//karta(:79)
2404         read (inp,'(a)') karta
2405         karta=ucase(karta)
2406       enddo
2407       card=card(:ilen(card)+1)//karta
2408       return
2409       end
2410 c----------------------------------------------------------------------------------
2411       subroutine readrst
2412       implicit real*8 (a-h,o-z)
2413       include 'DIMENSIONS'
2414       include 'COMMON.CHAIN'
2415       include 'COMMON.IOUNITS'
2416       include 'COMMON.MD'
2417       include 'COMMON.CONTROL'
2418       open(irest2,file=rest2name,status='unknown')
2419       read(irest2,*) totT,EK,potE,totE,t_bath
2420       do i=1,2*nres
2421          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2422       enddo
2423       do i=1,2*nres
2424          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2425       enddo
2426       if(usampl.or.homol_nset.gt.1) then
2427              read (irest2,*) iset
2428       endif
2429       close(irest2)
2430       return
2431       end
2432 c---------------------------------------------------------------------------------
2433       subroutine read_fragments
2434       implicit real*8 (a-h,o-z)
2435       include 'DIMENSIONS'
2436 #ifdef MPI
2437       include 'mpif.h'
2438 #endif
2439       include 'COMMON.SETUP'
2440       include 'COMMON.CHAIN'
2441       include 'COMMON.IOUNITS'
2442       include 'COMMON.MD'
2443       include 'COMMON.CONTROL'
2444       read(inp,*) nset,nfrag,npair,nfrag_back
2445       if(me.eq.king.or..not.out1file)
2446      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2447      &  " nfrag_back",nfrag_back
2448       do iset=1,nset
2449          read(inp,*) mset(iset)
2450        do i=1,nfrag
2451          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2452      &     qinfrag(i,iset)
2453          if(me.eq.king.or..not.out1file)
2454      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2455      &     ifrag(2,i,iset), qinfrag(i,iset)
2456        enddo
2457        do i=1,npair
2458         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2459      &    qinpair(i,iset)
2460         if(me.eq.king.or..not.out1file)
2461      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2462      &    ipair(2,i,iset), qinpair(i,iset)
2463        enddo
2464        do i=1,nfrag_back
2465         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2466      &     wfrag_back(3,i,iset),
2467      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2468         if(me.eq.king.or..not.out1file)
2469      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2470      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2471        enddo
2472       enddo
2473       return
2474       end
2475 c-------------------------------------------------------------------------------
2476       subroutine read_dist_constr
2477       implicit real*8 (a-h,o-z)
2478       include 'DIMENSIONS'
2479 #ifdef MPI
2480       include 'mpif.h'
2481 #endif
2482       include 'COMMON.SETUP'
2483       include 'COMMON.CONTROL'
2484       include 'COMMON.CHAIN'
2485       include 'COMMON.IOUNITS'
2486       include 'COMMON.SBRIDGE'
2487       integer ifrag_(2,100),ipair_(2,100)
2488       double precision wfrag_(100),wpair_(100)
2489       character*500 controlcard
2490 c      write (iout,*) "Calling read_dist_constr"
2491 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2492 c      call flush(iout)
2493       call card_concat(controlcard)
2494       call readi(controlcard,"NFRAG",nfrag_,0)
2495       call readi(controlcard,"NPAIR",npair_,0)
2496       call readi(controlcard,"NDIST",ndist_,0)
2497       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2498       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2499       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2500       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2501       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2502 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2503 c      write (iout,*) "IFRAG"
2504 c      do i=1,nfrag_
2505 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2506 c      enddo
2507 c      write (iout,*) "IPAIR"
2508 c      do i=1,npair_
2509 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2510 c      enddo
2511       if (.not.refstr .and. nfrag.gt.0) then
2512         write (iout,*)
2513      &  "ERROR: no reference structure to compute distance restraints"
2514         write (iout,*)
2515      &  "Restraints must be specified explicitly (NDIST=number)"
2516         stop
2517       endif
2518       if (nfrag.lt.2 .and. npair.gt.0) then
2519         write (iout,*) "ERROR: Less than 2 fragments specified",
2520      &   " but distance restraints between pairs requested"
2521         stop
2522       endif
2523       call flush(iout)
2524       do i=1,nfrag_
2525         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2526         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2527      &    ifrag_(2,i)=nstart_sup+nsup-1
2528 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2529         call flush(iout)
2530         if (wfrag_(i).gt.0.0d0) then
2531         do j=ifrag_(1,i),ifrag_(2,i)-1
2532           do k=j+1,ifrag_(2,i)
2533 c            write (iout,*) "j",j," k",k
2534             ddjk=dist(j,k)
2535             if (constr_dist.eq.1) then
2536             nhpb=nhpb+1
2537             ihpb(nhpb)=j
2538             jhpb(nhpb)=k
2539               dhpb(nhpb)=ddjk
2540             forcon(nhpb)=wfrag_(i)
2541             else if (constr_dist.eq.2) then
2542               if (ddjk.le.dist_cut) then
2543                 nhpb=nhpb+1
2544                 ihpb(nhpb)=j
2545                 jhpb(nhpb)=k
2546                 dhpb(nhpb)=ddjk
2547                 forcon(nhpb)=wfrag_(i)
2548               endif
2549             else
2550               nhpb=nhpb+1
2551               ihpb(nhpb)=j
2552               jhpb(nhpb)=k
2553               dhpb(nhpb)=ddjk
2554               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2555             endif
2556 #ifdef MPI
2557             if (.not.out1file .or. me.eq.king)
2558      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2559      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2560 #else
2561             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2562      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2563 #endif
2564           enddo
2565         enddo
2566         endif
2567       enddo
2568       do i=1,npair_
2569         if (wpair_(i).gt.0.0d0) then
2570         ii = ipair_(1,i)
2571         jj = ipair_(2,i)
2572         if (ii.gt.jj) then
2573           itemp=ii
2574           ii=jj
2575           jj=itemp
2576         endif
2577         do j=ifrag_(1,ii),ifrag_(2,ii)
2578           do k=ifrag_(1,jj),ifrag_(2,jj)
2579             nhpb=nhpb+1
2580             ihpb(nhpb)=j
2581             jhpb(nhpb)=k
2582             forcon(nhpb)=wpair_(i)
2583             dhpb(nhpb)=dist(j,k)
2584 #ifdef MPI
2585             if (.not.out1file .or. me.eq.king)
2586      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2587      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2588 #else
2589             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2590      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2591 #endif
2592           enddo
2593         enddo
2594         endif
2595       enddo
2596       do i=1,ndist_
2597         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2598      &     ibecarb(i),forcon(nhpb+1)
2599         if (forcon(nhpb+1).gt.0.0d0) then
2600           nhpb=nhpb+1
2601           if (ibecarb(i).gt.0) then
2602             ihpb(i)=ihpb(i)+nres
2603             jhpb(i)=jhpb(i)+nres
2604           endif
2605           if (dhpb(nhpb).eq.0.0d0)
2606      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2607         endif
2608       enddo
2609 #ifdef MPI
2610       if (.not.out1file .or. me.eq.king) then
2611 #endif
2612       do i=1,nhpb
2613           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2614      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2615       enddo
2616       call flush(iout)
2617 #ifdef MPI
2618       endif
2619 #endif
2620       return
2621       end
2622 c-------------------------------------------------------------------------------
2623
2624       subroutine read_constr_homology
2625
2626       include 'DIMENSIONS'
2627 #ifdef MPI
2628       include 'mpif.h'
2629 #endif
2630       include 'COMMON.SETUP'
2631       include 'COMMON.CONTROL'
2632       include 'COMMON.CHAIN'
2633       include 'COMMON.IOUNITS'
2634       include 'COMMON.MD'
2635       include 'COMMON.GEO'
2636       include 'COMMON.INTERACT'
2637       double precision odl_temp,sigma_odl_temp
2638       common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2639      &    sigma_odl_temp(maxres,maxres,max_template)
2640       character*2 kic2
2641       character*24 model_ki_dist, model_ki_angle
2642       character*500 controlcard
2643       character*3200 controlcard1
2644       integer ki, i, j, k, l
2645       logical lprn /.true./
2646
2647       call card_concat(controlcard)
2648       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2649       if (homol_nset.gt.1)then
2650          call card_concat(controlcard)
2651          read(controlcard,*) (waga_dist(i),i=1,homol_nset)
2652          call card_concat(controlcard)
2653          read(controlcard,*) (waga_angle(i),i=1,homol_nset)
2654          write(iout,*) "iset distance_weight angle_weight"
2655          do i=1,homol_nset
2656            write(iout,*) i,waga_dist(i),waga_angle(i)
2657          enddo
2658       else
2659        iset=1
2660        call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
2661        call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
2662       endif
2663
2664       write (iout,*) "nnt",nnt," nct",nct
2665       call flush(iout)
2666       lim_odl=0
2667       lim_dih=0
2668       do i=1,nres
2669         do j=i+2,nres
2670           do ki=1,constr_homology
2671             sigma_odl_temp(i,j,ki)=0.0d0
2672             odl_temp(i,j,ki)=0.0d0
2673           enddo
2674         enddo
2675       enddo
2676       do i=1,nres-3
2677         do ki=1,constr_homology
2678           dih(ki,i)=0.0d0
2679           sigma_dih(ki,i)=0.0d0
2680         enddo
2681       enddo
2682       do ki=1,constr_homology
2683           write(kic2,'(i2)') ki
2684           if (ki.le.9) kic2="0"//kic2(2:2)
2685
2686           model_ki_dist="model"//kic2//".dist"
2687           model_ki_angle="model"//kic2//".angle"
2688         open (ientin,file=model_ki_dist,status='old')
2689         do irec=1,maxdim !petla do czytania wiezow na odleglosc
2690           read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
2691      &       sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
2692           odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
2693           sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
2694      &     sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
2695         enddo
2696  1401 continue
2697         close (ientin)
2698         open (ientin,file=model_ki_angle,status='old')
2699         do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
2700           read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
2701      &      sigma_dih(ki,i+nnt-1)
2702           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
2703           sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
2704         enddo
2705  1402 continue
2706         close (ientin)
2707       enddo
2708       ii=0
2709       write (iout,*) "nnt",nnt," nct",nct
2710       do i=nnt,nct-2
2711         do j=i+2,nct
2712           ki=1
2713 c          write (iout,*) "i",i," j",j," constr_homology",constr_homology
2714           do while (ki.le.constr_homology .and.
2715      &        sigma_odl_temp(i,j,ki).le.0.0d0)
2716 c            write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
2717             ki=ki+1
2718           enddo
2719 c          write (iout,*) "ki",ki
2720           if (ki.gt.constr_homology) cycle
2721           ii=ii+1
2722           ires_homo(ii)=i
2723           jres_homo(ii)=j
2724           do ki=1,constr_homology
2725             odl(ki,ii)=odl_temp(i,j,ki)
2726             sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
2727           enddo
2728         enddo
2729       enddo
2730       lim_odl=ii
2731       if (constr_homology.gt.0) call homology_partition
2732 c Print restraints
2733       if (.not.lprn) return
2734       write (iout,*) "Distance restraints from templates"
2735       do ii=1,lim_odl
2736         write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
2737      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
2738       enddo
2739       write (iout,*) "Dihedral angle restraints from templates"
2740       do i=nnt,lim_dih
2741         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
2742      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2743       enddo
2744 c      write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
2745 c      write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
2746
2747
2748       return
2749       end
2750 c----------------------------------------------------------------------
2751
2752 #ifdef WINIFL
2753       subroutine flush(iu)
2754       return
2755       end
2756 #endif
2757 #ifdef AIX
2758       subroutine flush(iu)
2759       call flush_(iu)
2760       return
2761       end
2762 #endif
2763
2764 c------------------------------------------------------------------------------
2765       subroutine copy_to_tmp(source)
2766       include "DIMENSIONS"
2767       include "COMMON.IOUNITS"
2768       character*(*) source
2769       character* 256 tmpfile
2770       integer ilen
2771       external ilen
2772       logical ex
2773       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2774       inquire(file=tmpfile,exist=ex)
2775       if (ex) then
2776         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2777      &   " to temporary directory..."
2778         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2779         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2780       endif
2781       return
2782       end
2783 c------------------------------------------------------------------------------
2784       subroutine move_from_tmp(source)
2785       include "DIMENSIONS"
2786       include "COMMON.IOUNITS"
2787       character*(*) source
2788       integer ilen
2789       external ilen
2790       write (*,*) "Moving ",source(:ilen(source)),
2791      & " from temporary directory to working directory"
2792       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2793       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2794       return
2795       end
2796 c------------------------------------------------------------------------------
2797       subroutine random_init(seed)
2798 C
2799 C Initialize random number generator
2800 C
2801       implicit real*8 (a-h,o-z)
2802       include 'DIMENSIONS'
2803 #ifdef AMD64
2804       integer*8 iseedi8
2805 #endif
2806 #ifdef MPI
2807       include 'mpif.h'
2808       logical OKRandom, prng_restart
2809       real*8  r1
2810       integer iseed_array(4)
2811 #endif
2812       include 'COMMON.IOUNITS'
2813       include 'COMMON.TIME1'
2814       include 'COMMON.THREAD'
2815       include 'COMMON.SBRIDGE'
2816       include 'COMMON.CONTROL'
2817       include 'COMMON.MCM'
2818       include 'COMMON.MAP'
2819       include 'COMMON.HEADER'
2820 csa      include 'COMMON.CSA'
2821       include 'COMMON.CHAIN'
2822       include 'COMMON.MUCA'
2823       include 'COMMON.MD'
2824       include 'COMMON.FFIELD'
2825       include 'COMMON.SETUP'
2826       iseed=-dint(dabs(seed))
2827       if (iseed.eq.0) then
2828         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2829      &    'Random seed undefined. The program will stop.'
2830         write (*,'(/80(1h*)/20x,a/80(1h*))')
2831      &    'Random seed undefined. The program will stop.'
2832 #ifdef MPI
2833         call mpi_finalize(mpi_comm_world,ierr)
2834 #endif
2835         stop 'Bad random seed.'
2836       endif
2837 #ifdef MPI
2838       if (fg_rank.eq.0) then
2839       seed=seed*(me+1)+1
2840 #ifdef AMD64
2841       iseedi8=dint(seed)
2842       if(me.eq.king .or. .not. out1file)
2843      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2844       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2845       OKRandom = prng_restart(me,iseedi8)
2846 #else
2847       do i=1,4
2848        tmp=65536.0d0**(4-i)
2849        iseed_array(i) = dint(seed/tmp)
2850        seed=seed-iseed_array(i)*tmp
2851       enddo
2852       if(me.eq.king .or. .not. out1file)
2853      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2854      &                 (iseed_array(i),i=1,4)
2855       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2856      &                 (iseed_array(i),i=1,4)
2857       OKRandom = prng_restart(me,iseed_array)
2858 #endif
2859       if (OKRandom) then
2860         r1=ran_number(0.0D0,1.0D0)
2861         if(me.eq.king .or. .not. out1file)
2862      &   write (iout,*) 'ran_num',r1
2863         if (r1.lt.0.0d0) OKRandom=.false.
2864       endif
2865       if (.not.OKRandom) then
2866         write (iout,*) 'PRNG IS NOT WORKING!!!'
2867         print *,'PRNG IS NOT WORKING!!!'
2868         if (me.eq.0) then
2869          call flush(iout)
2870          call mpi_abort(mpi_comm_world,error_msg,ierr)
2871          stop
2872         else
2873          write (iout,*) 'too many processors for parallel prng'
2874          write (*,*) 'too many processors for parallel prng'
2875          call flush(iout)
2876          stop
2877         endif
2878       endif
2879       endif
2880 #else
2881       call vrndst(iseed)
2882       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2883 #endif
2884       return
2885       end