source/unres/src_CSA/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       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21 c      if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23 c      if (modecalc.eq.14) then
  24 c         call read_MDpar
  25 c         call read_REMDpar
  26 c      endif
  27 C Read MUCA control parameters if required
  28 c      if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39
  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 c      include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74 c      include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76 c      include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78 c      include 'COMMON.MUCA'
  79 c      include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 100       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 101       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 102       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 103       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 104       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 105       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 106       call reada(controlcard,'DRMS',drms,0.1D0)
 107       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 108        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 109        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 110        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 111        write (iout,'(a,f10.1)')'DRMS    = ',drms
 112        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 113        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 114       endif
 115       call readi(controlcard,'NZ_START',nz_start,0)
 116       call readi(controlcard,'NZ_END',nz_end,0)
 117       call readi(controlcard,'IZ_SC',iz_sc,0)
 118       timlim=60.0D0*timlim
 119       safety = 60.0d0*safety
 120       timem=timlim
 121       modecalc=0
 122       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 123       minim=(index(controlcard,'MINIMIZE').gt.0)
 124       dccart=(index(controlcard,'CART').gt.0)
 125       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 126       overlapsc=.not.overlapsc
 127       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 128       searchsc=.not.searchsc
 129       sideadd=(index(controlcard,'SIDEADD').gt.0)
 130       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 131       outpdb=(index(controlcard,'PDBOUT').gt.0)
 132       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 133       pdbref=(index(controlcard,'PDBREF').gt.0)
 134       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 135       indpdb=index(controlcard,'PDBSTART')
 136       extconf=(index(controlcard,'EXTCONF').gt.0)
 137       call readi(controlcard,'IPRINT',iprint,0)
 138       call readi(controlcard,'MAXGEN',maxgen,10000)
 139       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 140       call readi(controlcard,"KDIAG",kdiag,0)
 141       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
 142       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 143      & write (iout,*) "RESCALE_MODE",rescale_mode
 144       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 145       if (index(controlcard,'REGULAR').gt.0.0D0) then
 146         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 147         modecalc=1
 148         refstr=.true.
 149       endif
 150       if (index(controlcard,'CHECKGRAD').gt.0) then
 151         modecalc=5
 152         if (index(controlcard,'CART').gt.0) then
 153           icheckgrad=1
 154         elseif (index(controlcard,'CARINT').gt.0) then
 155           icheckgrad=2
 156         else
 157           icheckgrad=3
 158         endif
 159       elseif (index(controlcard,'THREAD').gt.0) then
 160         modecalc=2
 161         call readi(controlcard,'THREAD',nthread,0)
 162         if (nthread.gt.0) then
 163           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 164         else
 165           if (fg_rank.eq.0)
 166      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 167           stop 'Error termination in Read_Control.'
 168         endif
 169       else if (index(controlcard,'MCMA').gt.0) then
 170         modecalc=3
 171       else if (index(controlcard,'MCEE').gt.0) then
 172         modecalc=6
 173       else if (index(controlcard,'MULTCONF').gt.0) then
 174         modecalc=4
 175       else if (index(controlcard,'MAP').gt.0) then
 176         modecalc=7
 177         call readi(controlcard,'MAP',nmap,0)
 178       else if (index(controlcard,'CSA').gt.0) then
 179         modecalc=8
 180 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 181 crc
 182 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 183 crc
 184 crc        modecalc=9
 185 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 186 cfmc        modecalc=10
 187       else if (index(controlcard,'SOFTREG').gt.0) then
 188         modecalc=11
 189       else if (index(controlcard,'CHECK_BOND').gt.0) then
 190         modecalc=-1
 191       else if (index(controlcard,'TEST').gt.0) then
 192         modecalc=-2
 193       else if (index(controlcard,'MD').gt.0) then
 194         modecalc=12
 195       else if (index(controlcard,'RE ').gt.0) then
 196         modecalc=14
 197       endif
 198
 199       lmuca=index(controlcard,'MUCA').gt.0
 200       call readi(controlcard,'MUCADYN',mucadyn,0)
 201       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 202       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 203      & then
 204        write (iout,*) 'MUCADYN=',mucadyn
 205        write (iout,*) 'MUCASMOOTH=',muca_smooth
 206       endif
 207
 208       iscode=index(controlcard,'ONE_LETTER')
 209       indphi=index(controlcard,'PHI')
 210       indback=index(controlcard,'BACK')
 211       iranconf=index(controlcard,'RAND_CONF')
 212       i2ndstr=index(controlcard,'USE_SEC_PRED')
 213       gradout=index(controlcard,'GRADOUT').gt.0
 214       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 215
 216       if(me.eq.king.or..not.out1file)
 217      & write (iout,'(2a)') diagmeth(kdiag),
 218      &  ' routine used to diagonalize matrices.'
 219       return
 220       end
 221 c------------------------------------------------------------------------------
 222       subroutine molread
 223 C
 224 C Read molecular data.
 225 C
 226       implicit real*8 (a-h,o-z)
 227       include 'DIMENSIONS'
 228 #ifdef MPI
 229       include 'mpif.h'
 230       integer error_msg
 231 #endif
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.GEO'
 234       include 'COMMON.VAR'
 235       include 'COMMON.INTERACT'
 236       include 'COMMON.LOCAL'
 237       include 'COMMON.NAMES'
 238       include 'COMMON.CHAIN'
 239       include 'COMMON.FFIELD'
 240       include 'COMMON.SBRIDGE'
 241       include 'COMMON.HEADER'
 242       include 'COMMON.CONTROL'
 243 c      include 'COMMON.DBASE'
 244 c      include 'COMMON.THREAD'
 245       include 'COMMON.CONTACTS'
 246       include 'COMMON.TORCNSTR'
 247       include 'COMMON.TIME1'
 248       include 'COMMON.BOUNDS'
 249 c      include 'COMMON.MD'
 250 c      include 'COMMON.REMD'
 251       include 'COMMON.SETUP'
 252       character*4 sequence(maxres)
 253       integer rescode
 254       double precision x(maxvar)
 255       character*256 pdbfile
 256       character*320 weightcard
 257       character*80 weightcard_t,ucase
 258       dimension itype_pdb(maxres)
 259       common /pizda/ itype_pdb
 260       logical seq_comp,fail
 261       double precision energia(0:n_ene)
 262       integer ilen
 263       external ilen
 264 C
 265 C Body
 266 C
 267 C Read weights of the subsequent energy terms.
 268
 269        call card_concat(weightcard)
 270        call reada(weightcard,'WLONG',wlong,1.0D0)
 271        call reada(weightcard,'WSC',wsc,wlong)
 272        call reada(weightcard,'WSCP',wscp,wlong)
 273        call reada(weightcard,'WELEC',welec,1.0D0)
 274        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 275        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 276        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 277        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 278        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 279        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 280        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 281        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 282        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 283        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 284        call reada(weightcard,'WBOND',wbond,1.0D0)
 285        call reada(weightcard,'WTOR',wtor,1.0D0)
 286        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 287        call reada(weightcard,'WANG',wang,1.0D0)
 288        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 289 C     Juyong
 290        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 291        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 292        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 293        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 294 C
 295        call reada(weightcard,'SCAL14',scal14,0.4D0)
 296        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 297        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 298        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 299        call reada(weightcard,'TEMP0',temp0,300.0d0)
 300        if (index(weightcard,'SOFT').gt.0) ipot=6
 301 C 12/1/95 Added weight for the multi-body term WCORR
 302        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 303        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 304        weights(1)=wsc
 305        weights(2)=wscp
 306        weights(3)=welec
 307        weights(4)=wcorr
 308        weights(5)=wcorr5
 309        weights(6)=wcorr6
 310        weights(7)=wel_loc
 311        weights(8)=wturn3
 312        weights(9)=wturn4
 313        weights(10)=wturn6
 314        weights(11)=wang
 315        weights(12)=wscloc
 316        weights(13)=wtor
 317        weights(14)=wtor_d
 318        weights(15)=wstrain
 319        weights(16)=wvdwpp
 320        weights(17)=wbond
 321        weights(18)=scal14
 322        weights(21)=wsccor
 323 C     JUYONG
 324        weights(24)=wdfa_dist
 325        weights(25)=wdfa_tor
 326        weights(26)=wdfa_nei
 327        weights(27)=wdfa_beta
 328 C
 329
 330       if(me.eq.king.or..not.out1file)
 331      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 332      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 333      &  wturn4,wturn6,
 334      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 335
 336    10 format (/'Energy-term weights (unscaled):'//
 337      & 'WSCC=   ',f10.6,' (SC-SC)'/
 338      & 'WSCP=   ',f10.6,' (SC-p)'/
 339      & 'WELEC=  ',f10.6,' (p-p electr)'/
 340      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 341      & 'WBOND=  ',f10.6,' (stretching)'/
 342      & 'WANG=   ',f10.6,' (bending)'/
 343      & 'WSCLOC= ',f10.6,' (SC local)'/
 344      & 'WTOR=   ',f10.6,' (torsional)'/
 345      & 'WTORD=  ',f10.6,' (double torsional)'/
 346      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 347      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 348      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 349      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 350      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 351      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 352      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 353      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 354      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 355      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 356      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 357      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 358      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 359
 360       if(me.eq.king.or..not.out1file)then
 361        if (wcorr4.gt.0.0d0) then
 362         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 363      &   'between contact pairs of peptide groups'
 364         write (iout,'(2(a,f5.3/))')
 365      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 366      &  'Range of quenching the correlation terms:',2*delt_corr
 367        else if (wcorr.gt.0.0d0) then
 368         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 369      &   'between contact pairs of peptide groups'
 370        endif
 371        write (iout,'(a,f8.3)')
 372      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 373        write (iout,'(a,f8.3)')
 374      &  'General scaling factor of SC-p interactions:',scalscp
 375       endif
 376       r0_corr=cutoff_corr-delt_corr
 377       do i=1,20
 378         aad(i,1)=scalscp*aad(i,1)
 379         aad(i,2)=scalscp*aad(i,2)
 380         bad(i,1)=scalscp*bad(i,1)
 381         bad(i,2)=scalscp*bad(i,2)
 382       enddo
 383 c      call rescale_weights(t_bath)
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 386      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 387      &  wturn4,wturn6,
 388      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 389
 390    22 format (/'Energy-term weights (scaled):'//
 391      & 'WSCC=   ',f10.6,' (SC-SC)'/
 392      & 'WSCP=   ',f10.6,' (SC-p)'/
 393      & 'WELEC=  ',f10.6,' (p-p electr)'/
 394      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 395      & 'WBOND=  ',f10.6,' (stretching)'/
 396      & 'WANG=   ',f10.6,' (bending)'/
 397      & 'WSCLOC= ',f10.6,' (SC local)'/
 398      & 'WTOR=   ',f10.6,' (torsional)'/
 399      & 'WTORD=  ',f10.6,' (double torsional)'/
 400      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 401      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 402      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 403      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 404      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 405      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 406      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 407      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 408      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 409      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 410      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 411      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 412      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 413
 414       if(me.eq.king.or..not.out1file)
 415      & write (iout,*) "Reference temperature for weights calculation:",
 416      &  temp0
 417       call reada(weightcard,"D0CM",d0cm,3.78d0)
 418       call reada(weightcard,"AKCM",akcm,15.1d0)
 419       call reada(weightcard,"AKTH",akth,11.0d0)
 420       call reada(weightcard,"AKCT",akct,12.0d0)
 421       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 422       call reada(weightcard,"V2SS",v2ss,7.61d0)
 423       call reada(weightcard,"V3SS",v3ss,13.7d0)
 424       call reada(weightcard,"EBR",ebr,-5.50D0)
 425       if(me.eq.king.or..not.out1file) then
 426        write (iout,*) "Parameters of the SS-bond potential:"
 427        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 428      & " AKCT",akct
 429        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 430        write (iout,*) "EBR",ebr
 431        print *,'indpdb=',indpdb,' pdbref=',pdbref
 432       endif
 433       if (indpdb.gt.0 .or. pdbref) then
 434         read(inp,'(a)') pdbfile
 435         if(me.eq.king.or..not.out1file)
 436      &   write (iout,'(2a)') 'PDB data will be read from file ',
 437      &   pdbfile(:ilen(pdbfile))
 438         open(ipdbin,file=pdbfile,status='old',err=33)
 439         goto 34
 440   33    write (iout,'(a)') 'Error opening PDB file.'
 441         stop
 442   34    continue
 443 c        print *,'Begin reading pdb data'
 444         call readpdb
 445 c        print *,'Finished reading pdb data'
 446         if(me.eq.king.or..not.out1file)
 447      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 448      &   ' nstart_sup=',nstart_sup
 449         do i=1,nres
 450           itype_pdb(i)=itype(i)
 451         enddo
 452         close (ipdbin)
 453         nnt=nstart_sup
 454         nct=nstart_sup+nsup-1
 455         call contact(.false.,ncont_ref,icont_ref,co)
 456
 457         if (sideadd) then
 458          if(me.eq.king.or..not.out1file)
 459      &    write(iout,*)'Adding sidechains'
 460          maxsi=1000
 461          do i=2,nres-1
 462           iti=itype(i)
 463           if (iti.ne.10) then
 464             nsi=0
 465             fail=.true.
 466             do while (fail.and.nsi.le.maxsi)
 467 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 468               nsi=nsi+1
 469             enddo
 470             if(fail) write(iout,*)'Adding sidechain failed for res ',
 471      &              i,' after ',nsi,' trials'
 472           endif
 473          enddo
 474         endif
 475       endif
 476
 477       if (indpdb.eq.0) then
 478 C Read sequence if not taken from the pdb file.
 479         read (inp,*) nres
 480 c        print *,'nres=',nres
 481         if (iscode.gt.0) then
 482           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 483         else
 484           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 485         endif
 486 C Convert sequence to numeric code
 487         do i=1,nres
 488           itype(i)=rescode(i,sequence(i),iscode)
 489         enddo
 490 C Assign initial virtual bond lengths
 491         do i=2,nres
 492           vbld(i)=vbl
 493           vbld_inv(i)=vblinv
 494         enddo
 495         do i=2,nres-1
 496           vbld(i+nres)=dsc(itype(i))
 497           vbld_inv(i+nres)=dsc_inv(itype(i))
 498 c          write (iout,*) "i",i," itype",itype(i),
 499 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 500         enddo
 501       endif
 502 c      print *,nres
 503 c      print '(20i4)',(itype(i),i=1,nres)
 504       do i=1,nres
 505 #ifdef PROCOR
 506         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 507 #else
 508         if (itype(i).eq.21) then
 509 #endif
 510           itel(i)=0
 511 #ifdef PROCOR
 512         else if (itype(i+1).ne.20) then
 513 #else
 514         else if (itype(i).ne.20) then
 515 #endif
 516           itel(i)=1
 517         else
 518           itel(i)=2
 519         endif
 520       enddo
 521       if(me.eq.king.or..not.out1file)then
 522        write (iout,*) "ITEL"
 523        do i=1,nres-1
 524          write (iout,*) i,itype(i),itel(i)
 525        enddo
 526        print *,'Call Read_Bridge.'
 527       endif
 528       call read_bridge
 529 C 8/13/98 Set limits to generating the dihedral angles
 530       do i=1,nres
 531         phibound(1,i)=-pi
 532         phibound(2,i)=pi
 533       enddo
 534       read (inp,*) ndih_constr
 535       if (ndih_constr.gt.0) then
 536         read (inp,*) ftors
 537         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 538         if(me.eq.king.or..not.out1file)then
 539          write (iout,*)
 540      &   'There are',ndih_constr,' constraints on phi angles.'
 541          do i=1,ndih_constr
 542           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 543          enddo
 544         endif
 545         do i=1,ndih_constr
 546           phi0(i)=deg2rad*phi0(i)
 547           drange(i)=deg2rad*drange(i)
 548         enddo
 549         if(me.eq.king.or..not.out1file)
 550      &   write (iout,*) 'FTORS',ftors
 551         do i=1,ndih_constr
 552           ii = idih_constr(i)
 553           phibound(1,ii) = phi0(i)-drange(i)
 554           phibound(2,ii) = phi0(i)+drange(i)
 555         enddo
 556       endif
 557       nnt=1
 558 #ifdef MPI
 559       if (me.eq.king) then
 560 #endif
 561        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 562        do i=1,nres
 563          write (iout,'(a3,i5,2f10.1)')
 564      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 565        enddo
 566 #ifdef MP
 567       endif
 568 #endif
 569       nct=nres
 570 cd      print *,'NNT=',NNT,' NCT=',NCT
 571       if (itype(1).eq.21) nnt=2
 572       if (itype(nres).eq.21) nct=nct-1
 573
 574 C     Juyong:READ init_vars
 575 C     Initialize variables!
 576 C     Juyong:READ read_info
 577 C     READ fragment information!!
 578 C     both routines should be in dfa.F file!!
 579
 580       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0
 581      &      .and. wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 582        call init_dfa_vars
 583        print*, 'init_dfa_vars finished!'
 584        call read_dfa_info
 585        print*, 'read_dfa_info finished!'
 586       endif
 587 C
 588 C
 589
 590
 591       if (pdbref) then
 592         if(me.eq.king.or..not.out1file)
 593      &   write (iout,'(a,i3)') 'nsup=',nsup
 594         nstart_seq=nnt
 595         if (nsup.le.(nct-nnt+1)) then
 596           do i=0,nct-nnt+1-nsup
 597             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 598               nstart_seq=nnt+i
 599               goto 111
 600             endif
 601           enddo
 602           write (iout,'(a)')
 603      &            'Error - sequences to be superposed do not match.'
 604           stop
 605         else
 606           do i=0,nsup-(nct-nnt+1)
 607             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 608      &      then
 609               nstart_sup=nstart_sup+i
 610               nsup=nct-nnt+1
 611               goto 111
 612             endif
 613           enddo
 614           write (iout,'(a)')
 615      &            'Error - sequences to be superposed do not match.'
 616         endif
 617   111   continue
 618         if (nsup.eq.0) nsup=nct-nnt
 619         if (nstart_sup.eq.0) nstart_sup=nnt
 620         if (nstart_seq.eq.0) nstart_seq=nnt
 621         if(me.eq.king.or..not.out1file)
 622      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 623      &                 ' nstart_seq=',nstart_seq
 624       endif
 625 c--- Zscore rms -------
 626       if (nz_start.eq.0) nz_start=nnt
 627       if (nz_end.eq.0 .and. nsup.gt.0) then
 628         nz_end=nnt+nsup-1
 629       else if (nz_end.eq.0) then
 630         nz_end=nct
 631       endif
 632       if(me.eq.king.or..not.out1file)then
 633        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 634        write (iout,*) 'IZ_SC=',iz_sc
 635       endif
 636 c----------------------
 637       call init_int_table
 638       if (refstr) then
 639         if (.not.pdbref) then
 640           call read_angles(inp,*38)
 641           goto 39
 642    38     write (iout,'(a)') 'Error reading reference structure.'
 643 #ifdef MPI
 644           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 645           stop 'Error reading reference structure'
 646 #endif
 647    39     call chainbuild
 648           call setup_var
 649 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 650           nstart_sup=nnt
 651           nstart_seq=nnt
 652           nsup=nct-nnt+1
 653           do i=1,2*nres
 654             do j=1,3
 655               cref(j,i)=c(j,i)
 656             enddo
 657           enddo
 658           call contact(.true.,ncont_ref,icont_ref,co)
 659         endif
 660 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 661         call flush(iout)
 662         if (constr_dist.gt.0) call read_dist_constr
 663 c        write (iout,*) "After read_dist_constr nhpb",nhpb
 664         call hpb_partition
 665         if(me.eq.king.or..not.out1file)
 666      &   write (iout,*) 'Contact order:',co
 667         if (pdbref) then
 668         if(me.eq.king.or..not.out1file)
 669      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 670         do i=1,ncont_ref
 671           do j=1,2
 672             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 673           enddo
 674           if(me.eq.king.or..not.out1file)
 675      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 676      &     icont_ref(1,i),' ',
 677      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 678         enddo
 679         endif
 680       endif
 681       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 682      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 683      &    modecalc.ne.10) then
 684 C If input structure hasn't been supplied from the PDB file read or generate
 685 C initial geometry.
 686         if (iranconf.eq.0 .and. .not. extconf) then
 687           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 688      &     write (iout,'(a)') 'Initial geometry will be read in.'
 689           if (read_cart) then
 690             read(inp,'(8f10.5)',end=36,err=36)
 691      &       ((c(l,k),l=1,3),k=1,nres),
 692      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 693             call int_from_cart1(.false.)
 694             do i=1,nres-1
 695               do j=1,3
 696                 dc(j,i)=c(j,i+1)-c(j,i)
 697                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 698               enddo
 699             enddo
 700             do i=nnt,nct
 701               if (itype(i).ne.10) then
 702                 do j=1,3
 703                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 704                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 705                 enddo
 706               endif
 707             enddo
 708             return
 709           else
 710             call read_angles(inp,*36)
 711           endif
 712           goto 37
 713    36     write (iout,'(a)') 'Error reading angle file.'
 714 #ifdef MPI
 715           call mpi_finalize( MPI_COMM_WORLD,IERR )
 716 #endif
 717           stop 'Error reading angle file.'
 718    37     continue
 719         else if (extconf) then
 720          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 721      &    write (iout,'(a)') 'Extended chain initial geometry.'
 722          do i=3,nres
 723           theta(i)=90d0*deg2rad
 724          enddo
 725          do i=4,nres
 726           phi(i)=180d0*deg2rad
 727          enddo
 728          do i=2,nres-1
 729           alph(i)=110d0*deg2rad
 730          enddo
 731          do i=2,nres-1
 732           omeg(i)=-120d0*deg2rad
 733          enddo
 734         else
 735           if(me.eq.king.or..not.out1file)
 736      &     write (iout,'(a)') 'Random-generated initial geometry.'
 737
 738
 739 #ifdef MPI
 740           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 741      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 742 #endif
 743             do itrial=1,100
 744               itmp=1
 745 c              call gen_rand_conf(itmp,*30)
 746               goto 40
 747    30         write (iout,*) 'Failed to generate random conformation',
 748      &          ', itrial=',itrial
 749               write (*,*) 'Processor:',me,
 750      &          ' Failed to generate random conformation',
 751      &          ' itrial=',itrial
 752               call intout
 753
 754 #ifdef AIX
 755               call flush_(iout)
 756 #else
 757               call flush(iout)
 758 #endif
 759             enddo
 760             write (iout,'(a,i3,a)') 'Processor:',me,
 761      &        ' error in generating random conformation.'
 762             write (*,'(a,i3,a)') 'Processor:',me,
 763      &        ' error in generating random conformation.'
 764             call flush(iout)
 765 #ifdef MPI
 766             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 767    40       continue
 768           endif
 769 #else
 770           do itrial=1,100
 771             itmp=1
 772 c            call gen_rand_conf(itmp,*31)
 773             goto 40
 774    31       write (iout,*) 'Failed to generate random conformation',
 775      &        ', itrial=',itrial
 776             write (*,*) 'Failed to generate random conformation',
 777      &        ', itrial=',itrial
 778           enddo
 779           write (iout,'(a,i3,a)') 'Processor:',me,
 780      &      ' error in generating random conformation.'
 781           write (*,'(a,i3,a)') 'Processor:',me,
 782      &      ' error in generating random conformation.'
 783           stop
 784    40     continue
 785 #endif
 786         endif
 787       elseif (modecalc.eq.4) then
 788         read (inp,'(a)') intinname
 789         open (intin,file=intinname,status='old',err=333)
 790         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 791      &  write (iout,'(a)') 'intinname',intinname
 792         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 793         goto 334
 794   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 795 #ifdef MPI
 796         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 797 #endif
 798         stop 'Error opening angle file.'
 799   334   continue
 800
 801       endif
 802 C Generate distance constraints, if the PDB structure is to be regularized.
 803       if (nthread.gt.0) then
 804         call read_threadbase
 805       endif
 806       call setup_var
 807       if (me.eq.king .or. .not. out1file)
 808      & call intout
 809       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 810         write (iout,'(/a,i3,a)')
 811      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 812         write (iout,'(20i4)') (iss(i),i=1,ns)
 813         write (iout,'(/a/)') 'Pre-formed links are:'
 814         do i=1,nss
 815           i1=ihpb(i)-nres
 816           i2=jhpb(i)-nres
 817           it1=itype(i1)
 818           it2=itype(i2)
 819           if (me.eq.king.or..not.out1file)
 820      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 821      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 822      &    ebr,forcon(i)
 823         enddo
 824         write (iout,'(a)')
 825       endif
 826 c       if (i2ndstr.gt.0) call secstrp2dihc
 827 c      call geom_to_var(nvar,x)
 828 c      call etotal(energia(0))
 829 c      call enerprint(energia(0))
 830 c      call briefout(0,etot)
 831 c      stop
 832 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 833 cd    write (iout,'(a)') 'Variable list:'
 834 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 835 #ifdef MPI
 836       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 837      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 838      &  'Processor',myrank,': end reading molecular data.'
 839 #endif
 840       return
 841       end
 842 c--------------------------------------------------------------------------
 843       logical function seq_comp(itypea,itypeb,length)
 844       implicit none
 845       integer length,itypea(length),itypeb(length)
 846       integer i
 847       do i=1,length
 848         if (itypea(i).ne.itypeb(i)) then
 849           seq_comp=.false.
 850           return
 851         endif
 852       enddo
 853       seq_comp=.true.
 854       return
 855       end
 856 c-----------------------------------------------------------------------------
 857       subroutine read_bridge
 858 C Read information about disulfide bridges.
 859       implicit real*8 (a-h,o-z)
 860       include 'DIMENSIONS'
 861 #ifdef MPI
 862       include 'mpif.h'
 863 #endif
 864       include 'COMMON.IOUNITS'
 865       include 'COMMON.GEO'
 866       include 'COMMON.VAR'
 867       include 'COMMON.INTERACT'
 868       include 'COMMON.LOCAL'
 869       include 'COMMON.NAMES'
 870       include 'COMMON.CHAIN'
 871       include 'COMMON.FFIELD'
 872       include 'COMMON.SBRIDGE'
 873       include 'COMMON.HEADER'
 874       include 'COMMON.CONTROL'
 875 c      include 'COMMON.DBASE'
 876 c      include 'COMMON.THREAD'
 877       include 'COMMON.TIME1'
 878       include 'COMMON.SETUP'
 879 C Read bridging residues.
 880       read (inp,*) ns,(iss(i),i=1,ns)
 881       print *,'ns=',ns
 882       if(me.eq.king.or..not.out1file)
 883      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 884 C Check whether the specified bridging residues are cystines.
 885       do i=1,ns
 886         if (itype(iss(i)).ne.1) then
 887           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 888      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 889      &   ' can form a disulfide bridge?!!!'
 890           write (*,'(2a,i3,a)')
 891      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 892      &   ' can form a disulfide bridge?!!!'
 893 #ifdef MPI
 894          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 895          stop
 896 #endif
 897         endif
 898       enddo
 899 C Read preformed bridges.
 900       if (ns.gt.0) then
 901       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 902       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 903       if (nss.gt.0) then
 904         nhpb=nss
 905 C Check if the residues involved in bridges are in the specified list of
 906 C bridging residues.
 907         do i=1,nss
 908           do j=1,i-1
 909             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 910      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 911               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 912      &      ' contains residues present in other pairs.'
 913               write (*,'(a,i3,a)') 'Disulfide pair',i,
 914      &      ' contains residues present in other pairs.'
 915 #ifdef MPI
 916               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 917               stop
 918 #endif
 919             endif
 920           enddo
 921           do j=1,ns
 922             if (ihpb(i).eq.iss(j)) goto 10
 923           enddo
 924           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 925    10     continue
 926           do j=1,ns
 927             if (jhpb(i).eq.iss(j)) goto 20
 928           enddo
 929           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 930    20     continue
 931           dhpb(i)=dbr
 932           forcon(i)=fbr
 933         enddo
 934         do i=1,nss
 935           ihpb(i)=ihpb(i)+nres
 936           jhpb(i)=jhpb(i)+nres
 937         enddo
 938       endif
 939       endif
 940       return
 941       end
 942 c----------------------------------------------------------------------------
 943       subroutine read_x(kanal,*)
 944       implicit real*8 (a-h,o-z)
 945       include 'DIMENSIONS'
 946       include 'COMMON.GEO'
 947       include 'COMMON.VAR'
 948       include 'COMMON.CHAIN'
 949       include 'COMMON.IOUNITS'
 950       include 'COMMON.CONTROL'
 951       include 'COMMON.LOCAL'
 952       include 'COMMON.INTERACT'
 953 c Read coordinates from input
 954 c
 955       read(kanal,'(8f10.5)',end=10,err=10)
 956      &  ((c(l,k),l=1,3),k=1,nres),
 957      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 958       do j=1,3
 959         c(j,nres+1)=c(j,1)
 960         c(j,2*nres)=c(j,nres)
 961       enddo
 962       call int_from_cart1(.false.)
 963       do i=1,nres-1
 964         do j=1,3
 965           dc(j,i)=c(j,i+1)-c(j,i)
 966           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 967         enddo
 968       enddo
 969       do i=nnt,nct
 970         if (itype(i).ne.10) then
 971           do j=1,3
 972             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 973             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 974           enddo
 975         endif
 976       enddo
 977
 978       return
 979    10 return1
 980       end
 981 c------------------------------------------------------------------------------
 982       subroutine setup_var
 983       implicit real*8 (a-h,o-z)
 984       include 'DIMENSIONS'
 985       include 'COMMON.IOUNITS'
 986       include 'COMMON.GEO'
 987       include 'COMMON.VAR'
 988       include 'COMMON.INTERACT'
 989       include 'COMMON.LOCAL'
 990       include 'COMMON.NAMES'
 991       include 'COMMON.CHAIN'
 992       include 'COMMON.FFIELD'
 993       include 'COMMON.SBRIDGE'
 994       include 'COMMON.HEADER'
 995       include 'COMMON.CONTROL'
 996 c      include 'COMMON.DBASE'
 997 c      include 'COMMON.THREAD'
 998       include 'COMMON.TIME1'
 999 C Set up variable list.
1000       ntheta=nres-2
1001       nphi=nres-3
1002       nvar=ntheta+nphi
1003       nside=0
1004       do i=2,nres-1
1005         if (itype(i).ne.10) then
1006           nside=nside+1
1007           ialph(i,1)=nvar+nside
1008           ialph(nside,2)=i
1009         endif
1010       enddo
1011       if (indphi.gt.0) then
1012         nvar=nphi
1013       else if (indback.gt.0) then
1014         nvar=nphi+ntheta
1015       else
1016         nvar=nvar+2*nside
1017       endif
1018 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1019       return
1020       end
1021 c----------------------------------------------------------------------------
1022       subroutine gen_dist_constr
1023 C Generate CA distance constraints.
1024       implicit real*8 (a-h,o-z)
1025       include 'DIMENSIONS'
1026       include 'COMMON.IOUNITS'
1027       include 'COMMON.GEO'
1028       include 'COMMON.VAR'
1029       include 'COMMON.INTERACT'
1030       include 'COMMON.LOCAL'
1031       include 'COMMON.NAMES'
1032       include 'COMMON.CHAIN'
1033       include 'COMMON.FFIELD'
1034       include 'COMMON.SBRIDGE'
1035       include 'COMMON.HEADER'
1036       include 'COMMON.CONTROL'
1037 c      include 'COMMON.DBASE'
1038 c      include 'COMMON.THREAD'
1039       include 'COMMON.TIME1'
1040       dimension itype_pdb(maxres)
1041       common /pizda/ itype_pdb
1042       character*2 iden
1043 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1044 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1045 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1046 cd     & ' nsup',nsup
1047       do i=nstart_sup,nstart_sup+nsup-1
1048 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1049 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1050         do j=i+2,nstart_sup+nsup-1
1051           nhpb=nhpb+1
1052           ihpb(nhpb)=i+nstart_seq-nstart_sup
1053           jhpb(nhpb)=j+nstart_seq-nstart_sup
1054           forcon(nhpb)=weidis
1055           dhpb(nhpb)=dist(i,j)
1056         enddo
1057       enddo
1058 cd      write (iout,'(a)') 'Distance constraints:'
1059 cd      do i=nss+1,nhpb
1060 cd        ii=ihpb(i)
1061 cd        jj=jhpb(i)
1062 cd        iden='CA'
1063 cd        if (ii.gt.nres) then
1064 cd          iden='SC'
1065 cd          ii=ii-nres
1066 cd          jj=jj-nres
1067 cd        endif
1068 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1069 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1070 cd     &  dhpb(i),forcon(i)
1071 cd      enddo
1072       return
1073       end
1074 c----------------------------------------------------------------------------
1075       subroutine csaread
1076       implicit real*8 (a-h,o-z)
1077       include 'DIMENSIONS'
1078       include 'COMMON.IOUNITS'
1079       include 'COMMON.GEO'
1080       include 'COMMON.CSA'
1081       include 'COMMON.BANK'
1082       include 'COMMON.CONTROL'
1083       character*80 ucase
1084       character*620 mcmcard
1085       call card_concat(mcmcard)
1086
1087       call readi(mcmcard,'NCONF',nconf,50)
1088       call readi(mcmcard,'NADD',nadd,0)
1089       call readi(mcmcard,'JSTART',jstart,1)
1090       call readi(mcmcard,'JEND',jend,1)
1091       call readi(mcmcard,'NSTMAX',nstmax,500000)
1092       call readi(mcmcard,'N0',n0,1)
1093       call readi(mcmcard,'N1',n1,6)
1094       call readi(mcmcard,'N2',n2,4)
1095       call readi(mcmcard,'N3',n3,0)
1096       call readi(mcmcard,'N4',n4,0)
1097       call readi(mcmcard,'N5',n5,0)
1098       call readi(mcmcard,'N6',n6,10)
1099       call readi(mcmcard,'N7',n7,0)
1100       call readi(mcmcard,'N8',n8,0)
1101       call readi(mcmcard,'N9',n9,0)
1102       call readi(mcmcard,'N14',n14,0)
1103       call readi(mcmcard,'N15',n15,0)
1104       call readi(mcmcard,'N16',n16,0)
1105       call readi(mcmcard,'N17',n17,0)
1106       call readi(mcmcard,'N18',n18,0)
1107
1108       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1109
1110       call readi(mcmcard,'NDIFF',ndiff,2)
1111       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1112       call readi(mcmcard,'IS1',is1,1)
1113       call readi(mcmcard,'IS2',is2,8)
1114       call readi(mcmcard,'NRAN0',nran0,4)
1115       call readi(mcmcard,'NRAN1',nran1,2)
1116       call readi(mcmcard,'IRR',irr,1)
1117       call readi(mcmcard,'NSEED',nseed,20)
1118       call readi(mcmcard,'NTOTAL',ntotal,10000)
1119       call reada(mcmcard,'CUT1',cut1,2.0d0)
1120       call reada(mcmcard,'CUT2',cut2,5.0d0)
1121       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1122       call readi(mcmcard,'ICMAX',icmax,1)
1123       call readi(mcmcard,'NBANKM',nbankm,400)
1124       call readi(mcmcard,'IUCUT',iucut,2)
1125       call readi(mcmcard,'IRESTART',irestart,0)
1126 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1127       ntbankm=0
1128 c!bankt
1129       call reada(mcmcard,'DELE',dele,20.0d0)
1130       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1131       call readi(mcmcard,'IREF',iref,0)
1132       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1133       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1134       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1135       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1136       write (iout,*) "NCONF_IN",nconf_in
1137       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1138       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1139      &      " for torsional angles"
1140       return
1141       end
1142
1143       subroutine read_minim
1144       implicit real*8 (a-h,o-z)
1145       include 'DIMENSIONS'
1146       include 'COMMON.MINIM'
1147       include 'COMMON.IOUNITS'
1148       character*80 ucase
1149       character*320 minimcard
1150       call card_concat(minimcard)
1151       call readi(minimcard,'MAXMIN',maxmin,2000)
1152       call readi(minimcard,'MAXFUN',maxfun,5000)
1153       call readi(minimcard,'MINMIN',minmin,maxmin)
1154       call readi(minimcard,'MINFUN',minfun,maxmin)
1155       call reada(minimcard,'TOLF',tolf,1.0D-2)
1156       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1157       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1158       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1159       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1160       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1161      &         'Options in energy minimization:'
1162       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1163      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1164      & 'MinMin:',MinMin,' MinFun:',MinFun,
1165      & ' TolF:',TolF,' RTolF:',RTolF
1166       return
1167       end
1168 c----------------------------------------------------------------------------
1169       subroutine read_angles(kanal,*)
1170       implicit real*8 (a-h,o-z)
1171       include 'DIMENSIONS'
1172       include 'COMMON.GEO'
1173       include 'COMMON.VAR'
1174       include 'COMMON.CHAIN'
1175       include 'COMMON.IOUNITS'
1176       include 'COMMON.CONTROL'
1177 c Read angles from input
1178 c
1179        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1180        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1181        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1182        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1183
1184        do i=1,nres
1185 c 9/7/01 avoid 180 deg valence angle
1186         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1187 c
1188         theta(i)=deg2rad*theta(i)
1189         phi(i)=deg2rad*phi(i)
1190         alph(i)=deg2rad*alph(i)
1191         omeg(i)=deg2rad*omeg(i)
1192        enddo
1193       return
1194    10 return1
1195       end
1196 c----------------------------------------------------------------------------
1197       subroutine reada(rekord,lancuch,wartosc,default)
1198       implicit none
1199       character*(*) rekord,lancuch
1200       double precision wartosc,default
1201       integer ilen,iread
1202       external ilen
1203       iread=index(rekord,lancuch)
1204       if (iread.eq.0) then
1205         wartosc=default
1206         return
1207       endif
1208       iread=iread+ilen(lancuch)+1
1209       read (rekord(iread:),*,err=10,end=10) wartosc
1210       return
1211   10  wartosc=default
1212       return
1213       end
1214 c----------------------------------------------------------------------------
1215       subroutine readi(rekord,lancuch,wartosc,default)
1216       implicit none
1217       character*(*) rekord,lancuch
1218       integer wartosc,default
1219       integer ilen,iread
1220       external ilen
1221       iread=index(rekord,lancuch)
1222       if (iread.eq.0) then
1223         wartosc=default
1224         return
1225       endif
1226       iread=iread+ilen(lancuch)+1
1227       read (rekord(iread:),*,err=10,end=10) wartosc
1228       return
1229   10  wartosc=default
1230       return
1231       end
1232 c----------------------------------------------------------------------------
1233       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1234       implicit none
1235       integer dim,i
1236       integer tablica(dim),default
1237       character*(*) rekord,lancuch
1238       character*80 aux
1239       integer ilen,iread
1240       external ilen
1241       do i=1,dim
1242         tablica(i)=default
1243       enddo
1244       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1245       if (iread.eq.0) return
1246       iread=iread+ilen(lancuch)+1
1247       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1248    10 return
1249       end
1250 c----------------------------------------------------------------------------
1251       subroutine multreada(rekord,lancuch,tablica,dim,default)
1252       implicit none
1253       integer dim,i
1254       double precision tablica(dim),default
1255       character*(*) rekord,lancuch
1256       character*80 aux
1257       integer ilen,iread
1258       external ilen
1259       do i=1,dim
1260         tablica(i)=default
1261       enddo
1262       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1263       if (iread.eq.0) return
1264       iread=iread+ilen(lancuch)+1
1265       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1266    10 return
1267       end
1268 c----------------------------------------------------------------------------
1269       subroutine openunits
1270       implicit real*8 (a-h,o-z)
1271       include 'DIMENSIONS'
1272 #ifdef MPI
1273       include 'mpif.h'
1274       character*16 form,nodename
1275       integer nodelen
1276 #endif
1277       include 'COMMON.SETUP'
1278       include 'COMMON.IOUNITS'
1279 c      include 'COMMON.MD'
1280       include 'COMMON.CONTROL'
1281       integer lenpre,lenpot,ilen,lentmp
1282       external ilen
1283       character*3 out1file_text,ucase
1284       character*3 ll
1285       external ucase
1286 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1287       call getenv_loc("PREFIX",prefix)
1288       pref_orig = prefix
1289       call getenv_loc("POT",pot)
1290       call getenv_loc("DIRTMP",tmpdir)
1291       call getenv_loc("CURDIR",curdir)
1292       call getenv_loc("OUT1FILE",out1file_text)
1293 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1294       out1file_text=ucase(out1file_text)
1295       if (out1file_text(1:1).eq."Y") then
1296         out1file=.true.
1297       else
1298         out1file=fg_rank.gt.0
1299       endif
1300       lenpre=ilen(prefix)
1301       lenpot=ilen(pot)
1302       lentmp=ilen(tmpdir)
1303       if (lentmp.gt.0) then
1304           write (*,'(80(1h!))')
1305           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1306           write (*,'(80(1h!))')
1307           write (*,*)"All output files will be on node /tmp directory."
1308 #ifdef MPI
1309         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1310         if (me.eq.king) then
1311           write (*,*) "The master node is ",nodename
1312         else if (fg_rank.eq.0) then
1313           write (*,*) "I am the CG slave node ",nodename
1314         else
1315           write (*,*) "I am the FG slave node ",nodename
1316         endif
1317 #endif
1318         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1319         lenpre = lentmp+lenpre+1
1320       endif
1321       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1322 C Get the names and open the input files
1323 #if defined(WINIFL) || defined(WINPGI)
1324       open(1,file=pref_orig(:ilen(pref_orig))//
1325      &  '.inp',status='old',readonly,shared)
1326        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1327 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1328 C Get parameter filenames and open the parameter files.
1329       call getenv_loc('BONDPAR',bondname)
1330       open (ibond,file=bondname,status='old',readonly,shared)
1331       call getenv_loc('THETPAR',thetname)
1332       open (ithep,file=thetname,status='old',readonly,shared)
1333 #ifndef CRYST_THETA
1334       call getenv_loc('THETPARPDB',thetname_pdb)
1335       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1336 #endif
1337       call getenv_loc('ROTPAR',rotname)
1338       open (irotam,file=rotname,status='old',readonly,shared)
1339 #ifndef CRYST_SC
1340       call getenv_loc('ROTPARPDB',rotname_pdb)
1341       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1342 #endif
1343       call getenv_loc('TORPAR',torname)
1344       open (itorp,file=torname,status='old',readonly,shared)
1345       call getenv_loc('TORDPAR',tordname)
1346       open (itordp,file=tordname,status='old',readonly,shared)
1347       call getenv_loc('FOURIER',fouriername)
1348       open (ifourier,file=fouriername,status='old',readonly,shared)
1349       call getenv_loc('ELEPAR',elename)
1350       open (ielep,file=elename,status='old',readonly,shared)
1351       call getenv_loc('SIDEPAR',sidename)
1352       open (isidep,file=sidename,status='old',readonly,shared)
1353 #elif (defined CRAY) || (defined AIX)
1354       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1355      &  action='read')
1356 c      print *,"Processor",myrank," opened file 1"
1357       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1358 c      print *,"Processor",myrank," opened file 9"
1359 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1360 C Get parameter filenames and open the parameter files.
1361       call getenv_loc('BONDPAR',bondname)
1362       open (ibond,file=bondname,status='old',action='read')
1363 c      print *,"Processor",myrank," opened file IBOND"
1364       call getenv_loc('THETPAR',thetname)
1365       open (ithep,file=thetname,status='old',action='read')
1366 c      print *,"Processor",myrank," opened file ITHEP"
1367 #ifndef CRYST_THETA
1368       call getenv_loc('THETPARPDB',thetname_pdb)
1369       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1370 #endif
1371       call getenv_loc('ROTPAR',rotname)
1372       open (irotam,file=rotname,status='old',action='read')
1373 c      print *,"Processor",myrank," opened file IROTAM"
1374 #ifndef CRYST_SC
1375       call getenv_loc('ROTPARPDB',rotname_pdb)
1376       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1377 #endif
1378       call getenv_loc('TORPAR',torname)
1379       open (itorp,file=torname,status='old',action='read')
1380 c      print *,"Processor",myrank," opened file ITORP"
1381       call getenv_loc('TORDPAR',tordname)
1382       open (itordp,file=tordname,status='old',action='read')
1383 c      print *,"Processor",myrank," opened file ITORDP"
1384       call getenv_loc('SCCORPAR',sccorname)
1385       open (isccor,file=sccorname,status='old',action='read')
1386 c      print *,"Processor",myrank," opened file ISCCOR"
1387       call getenv_loc('FOURIER',fouriername)
1388       open (ifourier,file=fouriername,status='old',action='read')
1389 c      print *,"Processor",myrank," opened file IFOURIER"
1390       call getenv_loc('ELEPAR',elename)
1391       open (ielep,file=elename,status='old',action='read')
1392 c      print *,"Processor",myrank," opened file IELEP"
1393       call getenv_loc('SIDEPAR',sidename)
1394       open (isidep,file=sidename,status='old',action='read')
1395 c      print *,"Processor",myrank," opened file ISIDEP"
1396 c      print *,"Processor",myrank," opened parameter files"
1397 #elif (defined G77)
1398       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1399       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1400 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1401 C Get parameter filenames and open the parameter files.
1402       call getenv_loc('BONDPAR',bondname)
1403       open (ibond,file=bondname,status='old')
1404       call getenv_loc('THETPAR',thetname)
1405       open (ithep,file=thetname,status='old')
1406 #ifndef CRYST_THETA
1407       call getenv_loc('THETPARPDB',thetname_pdb)
1408       open (ithep_pdb,file=thetname_pdb,status='old')
1409 #endif
1410       call getenv_loc('ROTPAR',rotname)
1411       open (irotam,file=rotname,status='old')
1412 #ifndef CRYST_SC
1413       call getenv_loc('ROTPARPDB',rotname_pdb)
1414       open (irotam_pdb,file=rotname_pdb,status='old')
1415 #endif
1416       call getenv_loc('TORPAR',torname)
1417       open (itorp,file=torname,status='old')
1418       call getenv_loc('TORDPAR',tordname)
1419       open (itordp,file=tordname,status='old')
1420       call getenv_loc('SCCORPAR',sccorname)
1421       open (isccor,file=sccorname,status='old')
1422       call getenv_loc('FOURIER',fouriername)
1423       open (ifourier,file=fouriername,status='old')
1424       call getenv_loc('ELEPAR',elename)
1425       open (ielep,file=elename,status='old')
1426       call getenv_loc('SIDEPAR',sidename)
1427       open (isidep,file=sidename,status='old')
1428 #else
1429       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1430      &  readonly)
1431        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1432 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1433 C Get parameter filenames and open the parameter files.
1434       call getenv_loc('BONDPAR',bondname)
1435       open (ibond,file=bondname,status='old',readonly)
1436       call getenv_loc('THETPAR',thetname)
1437       open (ithep,file=thetname,status='old',readonly)
1438 #ifndef CRYST_THETA
1439       call getenv_loc('THETPARPDB',thetname_pdb)
1440       print *,"thetname_pdb ",thetname_pdb
1441       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1442       print *,ithep_pdb," opened"
1443 #endif
1444       call getenv_loc('ROTPAR',rotname)
1445       open (irotam,file=rotname,status='old',readonly)
1446 #ifndef CRYST_SC
1447       call getenv_loc('ROTPARPDB',rotname_pdb)
1448       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1449 #endif
1450       call getenv_loc('TORPAR',torname)
1451       open (itorp,file=torname,status='old',readonly)
1452       call getenv_loc('TORDPAR',tordname)
1453       open (itordp,file=tordname,status='old',readonly)
1454       call getenv_loc('SCCORPAR',sccorname)
1455       open (isccor,file=sccorname,status='old',readonly)
1456       call getenv_loc('FOURIER',fouriername)
1457       open (ifourier,file=fouriername,status='old',readonly)
1458       call getenv_loc('ELEPAR',elename)
1459       open (ielep,file=elename,status='old',readonly)
1460       call getenv_loc('SIDEPAR',sidename)
1461       open (isidep,file=sidename,status='old',readonly)
1462 #endif
1463 #ifndef OLDSCP
1464 C
1465 C 8/9/01 In the newest version SCp interaction constants are read from a file
1466 C Use -DOLDSCP to use hard-coded constants instead.
1467 C
1468       call getenv_loc('SCPPAR',scpname)
1469 #if defined(WINIFL) || defined(WINPGI)
1470       open (iscpp,file=scpname,status='old',readonly,shared)
1471 #elif (defined CRAY)  || (defined AIX)
1472       open (iscpp,file=scpname,status='old',action='read')
1473 #elif (defined G77)
1474       open (iscpp,file=scpname,status='old')
1475 #else
1476       open (iscpp,file=scpname,status='old',readonly)
1477 #endif
1478 #endif
1479       call getenv_loc('PATTERN',patname)
1480 #if defined(WINIFL) || defined(WINPGI)
1481       open (icbase,file=patname,status='old',readonly,shared)
1482 #elif (defined CRAY)  || (defined AIX)
1483       open (icbase,file=patname,status='old',action='read')
1484 #elif (defined G77)
1485       open (icbase,file=patname,status='old')
1486 #else
1487       open (icbase,file=patname,status='old',readonly)
1488 #endif
1489 #ifdef MPI
1490 C Open output file only for CG processes
1491 c      print *,"Processor",myrank," fg_rank",fg_rank
1492       if (fg_rank.eq.0) then
1493
1494       if (nodes.eq.1) then
1495         npos=3
1496       else
1497         npos = dlog10(dfloat(nodes-1))+1
1498       endif
1499       if (npos.lt.3) npos=3
1500       write (liczba,'(i1)') npos
1501       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1502      &  //')'
1503       write (liczba,form) me
1504       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1505      &  liczba(:ilen(liczba))
1506       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1507      &  //'.int'
1508       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1509      &  //'.pdb'
1510       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1511      &  liczba(:ilen(liczba))//'.mol2'
1512       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1513      &  liczba(:ilen(liczba))//'.stat'
1514       if (lentmp.gt.0)
1515      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1516      &      //liczba(:ilen(liczba))//'.stat')
1517       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1518      &  //'.rst'
1519 c      if(usampl) then
1520 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1521 c     & liczba(:ilen(liczba))//'.const'
1522 c      endif
1523
1524       endif
1525 #else
1526       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1527       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1528       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1529       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1530       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1531       if (lentmp.gt.0)
1532      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1533      &    '.stat')
1534       rest2name=prefix(:ilen(prefix))//'.rst'
1535 c      if(usampl) then
1536 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1537 c      endif
1538 #endif
1539 #if defined(AIX) || defined(PGI)
1540       if (me.eq.king .or. .not. out1file)
1541      &   open(iout,file=outname,status='unknown')
1542 #ifdef DEBUG
1543       if (fg_rank.gt.0) then
1544         write (liczba,'(i3.3)') myrank/nfgtasks
1545         write (ll,'(bz,i3.3)') fg_rank
1546         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1547      &   status='unknown')
1548       endif
1549 #endif
1550       if(me.eq.king) then
1551        open(igeom,file=intname,status='unknown',position='append')
1552        open(ipdb,file=pdbname,status='unknown')
1553        open(imol2,file=mol2name,status='unknown')
1554        open(istat,file=statname,status='unknown',position='append')
1555       else
1556 c1out       open(iout,file=outname,status='unknown')
1557       endif
1558 #else
1559       if (me.eq.king .or. .not.out1file)
1560      &    open(iout,file=outname,status='unknown')
1561 #ifdef DEBUG
1562       if (fg_rank.gt.0) then
1563         write (liczba,'(i3.3)') myrank/nfgtasks
1564         write (ll,'(bz,i3.3)') fg_rank
1565         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1566      &   status='unknown')
1567       endif
1568 #endif
1569       if(me.eq.king) then
1570        open(igeom,file=intname,status='unknown',access='append')
1571        open(ipdb,file=pdbname,status='unknown')
1572        open(imol2,file=mol2name,status='unknown')
1573        open(istat,file=statname,status='unknown',access='append')
1574       else
1575 c1out       open(iout,file=outname,status='unknown')
1576       endif
1577 #endif
1578       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1579       csa_seed=prefix(:lenpre)//'.CSA.seed'
1580       csa_history=prefix(:lenpre)//'.CSA.history'
1581       csa_bank=prefix(:lenpre)//'.CSA.bank'
1582       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1583       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1584       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1585 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1586       csa_int=prefix(:lenpre)//'.int'
1587       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1588       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1589       csa_in=prefix(:lenpre)//'.CSA.in'
1590 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1591 C Write file names
1592       if (me.eq.king)then
1593       write (iout,'(80(1h-))')
1594       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1595       write (iout,'(80(1h-))')
1596       write (iout,*) "Input file                      : ",
1597      &  pref_orig(:ilen(pref_orig))//'.inp'
1598       write (iout,*) "Output file                     : ",
1599      &  outname(:ilen(outname))
1600       write (iout,*)
1601       write (iout,*) "Sidechain potential file        : ",
1602      &  sidename(:ilen(sidename))
1603 #ifndef OLDSCP
1604       write (iout,*) "SCp potential file              : ",
1605      &  scpname(:ilen(scpname))
1606 #endif
1607       write (iout,*) "Electrostatic potential file    : ",
1608      &  elename(:ilen(elename))
1609       write (iout,*) "Cumulant coefficient file       : ",
1610      &  fouriername(:ilen(fouriername))
1611       write (iout,*) "Torsional parameter file        : ",
1612      &  torname(:ilen(torname))
1613       write (iout,*) "Double torsional parameter file : ",
1614      &  tordname(:ilen(tordname))
1615       write (iout,*) "SCCOR parameter file : ",
1616      &  sccorname(:ilen(sccorname))
1617       write (iout,*) "Bond & inertia constant file    : ",
1618      &  bondname(:ilen(bondname))
1619       write (iout,*) "Bending parameter file          : ",
1620      &  thetname(:ilen(thetname))
1621       write (iout,*) "Rotamer parameter file          : ",
1622      &  rotname(:ilen(rotname))
1623       write (iout,*) "Threading database              : ",
1624      &  patname(:ilen(patname))
1625       if (lentmp.ne.0)
1626      &write (iout,*)" DIRTMP                          : ",
1627      &  tmpdir(:lentmp)
1628       write (iout,'(80(1h-))')
1629       endif
1630       return
1631       end
1632 c----------------------------------------------------------------------------
1633       subroutine card_concat(card)
1634       implicit real*8 (a-h,o-z)
1635       include 'DIMENSIONS'
1636       include 'COMMON.IOUNITS'
1637       character*(*) card
1638       character*80 karta,ucase
1639       external ilen
1640       read (inp,'(a)') karta
1641       karta=ucase(karta)
1642       card=' '
1643       do while (karta(80:80).eq.'&')
1644         card=card(:ilen(card)+1)//karta(:79)
1645         read (inp,'(a)') karta
1646         karta=ucase(karta)
1647       enddo
1648       card=card(:ilen(card)+1)//karta
1649       return
1650       end
1651
1652       subroutine read_dist_constr
1653       implicit real*8 (a-h,o-z)
1654       include 'DIMENSIONS'
1655 #ifdef MPI
1656       include 'mpif.h'
1657 #endif
1658       include 'COMMON.SETUP'
1659       include 'COMMON.CONTROL'
1660       include 'COMMON.CHAIN'
1661       include 'COMMON.IOUNITS'
1662       include 'COMMON.SBRIDGE'
1663       integer ifrag_(2,100),ipair_(2,100)
1664       double precision wfrag_(100),wpair_(100)
1665       character*500 controlcard
1666       write (iout,*) "Calling read_dist_constr"
1667       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1668       call flush(iout)
1669       call card_concat(controlcard)
1670       call readi(controlcard,"NFRAG",nfrag_,0)
1671       call readi(controlcard,"NPAIR",npair_,0)
1672       call readi(controlcard,"NDIST",ndist_,0)
1673       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1674       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1675       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1676       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1677       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1678 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1679 c      write (iout,*) "IFRAG"
1680 c      do i=1,nfrag_
1681 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1682 c      enddo
1683 c      write (iout,*) "IPAIR"
1684 c      do i=1,npair_
1685 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1686 c      enddo
1687       call flush(iout)
1688       do i=1,nfrag_
1689         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1690         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1691      &    ifrag_(2,i)=nstart_sup+nsup-1
1692 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1693         call flush(iout)
1694         if (wfrag_(i).gt.0.0d0) then
1695         do j=ifrag_(1,i),ifrag_(2,i)-1
1696           do k=j+1,ifrag_(2,i)
1697             write (iout,*) "j",j," k",k
1698             ddjk=dist(j,k)
1699             if (constr_dist.eq.1) then
1700             nhpb=nhpb+1
1701             ihpb(nhpb)=j
1702             jhpb(nhpb)=k
1703               dhpb(nhpb)=ddjk
1704             forcon(nhpb)=wfrag_(i)
1705             else if (constr_dist.eq.2) then
1706               if (ddjk.le.dist_cut) then
1707                 nhpb=nhpb+1
1708                 ihpb(nhpb)=j
1709                 jhpb(nhpb)=k
1710                 dhpb(nhpb)=ddjk
1711                 forcon(nhpb)=wfrag_(i)
1712               endif
1713             else
1714               nhpb=nhpb+1
1715               ihpb(nhpb)=j
1716               jhpb(nhpb)=k
1717               dhpb(nhpb)=ddjk
1718               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1719             endif
1720 #ifdef MPI
1721             if (.not.out1file .or. me.eq.king)
1722      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1723      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1724 #else
1725             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1726      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1727 #endif
1728           enddo
1729         enddo
1730         endif
1731       enddo
1732       do i=1,npair_
1733         if (wpair_(i).gt.0.0d0) then
1734         ii = ipair_(1,i)
1735         jj = ipair_(2,i)
1736         if (ii.gt.jj) then
1737           itemp=ii
1738           ii=jj
1739           jj=itemp
1740         endif
1741         do j=ifrag_(1,ii),ifrag_(2,ii)
1742           do k=ifrag_(1,jj),ifrag_(2,jj)
1743             nhpb=nhpb+1
1744             ihpb(nhpb)=j
1745             jhpb(nhpb)=k
1746             forcon(nhpb)=wpair_(i)
1747             dhpb(nhpb)=dist(j,k)
1748 #ifdef MPI
1749             if (.not.out1file .or. me.eq.king)
1750      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1751      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1752 #else
1753             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1754      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1755 #endif
1756           enddo
1757         enddo
1758         endif
1759       enddo
1760       do i=1,ndist_
1761         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1762         if (forcon(nhpb+1).gt.0.0d0) then
1763           nhpb=nhpb+1
1764           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1765 #ifdef MPI
1766           if (.not.out1file .or. me.eq.king)
1767      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1768      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1769 #else
1770           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1771      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1772 #endif
1773         endif
1774       enddo
1775       call flush(iout)
1776       return
1777       end
1778 c-------------------------------------------------------------------------------
1779 #ifdef WINIFL
1780       subroutine flush(iu)
1781       return
1782       end
1783 #endif
1784 #ifdef AIX
1785       subroutine flush(iu)
1786       call flush_(iu)
1787       return
1788       end
1789 #endif
1790 c------------------------------------------------------------------------------
1791       subroutine copy_to_tmp(source)
1792       include "DIMENSIONS"
1793       include "COMMON.IOUNITS"
1794       character*(*) source
1795       character* 256 tmpfile
1796       integer ilen
1797       external ilen
1798       logical ex
1799       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1800       inquire(file=tmpfile,exist=ex)
1801       if (ex) then
1802         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1803      &   " to temporary directory..."
1804         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1805         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1806       endif
1807       return
1808       end
1809 c------------------------------------------------------------------------------
1810       subroutine move_from_tmp(source)
1811       include "DIMENSIONS"
1812       include "COMMON.IOUNITS"
1813       character*(*) source
1814       integer ilen
1815       external ilen
1816       write (*,*) "Moving ",source(:ilen(source)),
1817      & " from temporary directory to working directory"
1818       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1819       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1820       return
1821       end
1822 c------------------------------------------------------------------------------
1823       subroutine random_init(seed)
1824 C
1825 C Initialize random number generator
1826 C
1827       implicit real*8 (a-h,o-z)
1828       include 'DIMENSIONS'
1829 #ifdef AMD64
1830       integer*8 iseedi8
1831 #endif
1832 #ifdef MPI
1833       include 'mpif.h'
1834       logical OKRandom, prng_restart
1835       real*8  r1
1836       integer iseed_array(4)
1837 #endif
1838       include 'COMMON.IOUNITS'
1839       include 'COMMON.TIME1'
1840 c      include 'COMMON.THREAD'
1841       include 'COMMON.SBRIDGE'
1842       include 'COMMON.CONTROL'
1843       include 'COMMON.MCM'
1844 c      include 'COMMON.MAP'
1845       include 'COMMON.HEADER'
1846 c      include 'COMMON.CSA'
1847       include 'COMMON.CHAIN'
1848 c      include 'COMMON.MUCA'
1849 c      include 'COMMON.MD'
1850       include 'COMMON.FFIELD'
1851       include 'COMMON.SETUP'
1852       iseed=-dint(dabs(seed))
1853       if (iseed.eq.0) then
1854         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1855      &    'Random seed undefined. The program will stop.'
1856         write (*,'(/80(1h*)/20x,a/80(1h*))')
1857      &    'Random seed undefined. The program will stop.'
1858 #ifdef MPI
1859         call mpi_finalize(mpi_comm_world,ierr)
1860 #endif
1861         stop 'Bad random seed.'
1862       endif
1863 #ifdef MPI
1864       if (fg_rank.eq.0) then
1865       seed=seed*(me+1)+1
1866 #ifdef AMD64
1867       iseedi8=dint(seed)
1868       if(me.eq.king .or. .not. out1file)
1869      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1870       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1871       OKRandom = prng_restart(me,iseedi8)
1872 #else
1873       do i=1,4
1874        tmp=65536.0d0**(4-i)
1875        iseed_array(i) = dint(seed/tmp)
1876        seed=seed-iseed_array(i)*tmp
1877       enddo
1878       if(me.eq.king .or. .not. out1file)
1879      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1880      &                 (iseed_array(i),i=1,4)
1881       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1882      &                 (iseed_array(i),i=1,4)
1883       OKRandom = prng_restart(me,iseed_array)
1884 #endif
1885       if (OKRandom) then
1886         r1=ran_number(0.0D0,1.0D0)
1887         if(me.eq.king .or. .not. out1file)
1888      &   write (iout,*) 'ran_num',r1
1889         if (r1.lt.0.0d0) OKRandom=.false.
1890       endif
1891       if (.not.OKRandom) then
1892         write (iout,*) 'PRNG IS NOT WORKING!!!'
1893         print *,'PRNG IS NOT WORKING!!!'
1894         if (me.eq.0) then
1895          call flush(iout)
1896          call mpi_abort(mpi_comm_world,error_msg,ierr)
1897          stop
1898         else
1899          write (iout,*) 'too many processors for parallel prng'
1900          write (*,*) 'too many processors for parallel prng'
1901          call flush(iout)
1902          stop
1903         endif
1904       endif
1905       endif
1906 #else
1907       call vrndst(iseed)
1908       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1909 #endif
1910       return
1911       end