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