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