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