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,1)
 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       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       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1158      &         'Options in energy minimization:'
1159       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1160      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1161      & 'MinMin:',MinMin,' MinFun:',MinFun,
1162      & ' TolF:',TolF,' RTolF:',RTolF
1163       return
1164       end
1165 c----------------------------------------------------------------------------
1166       subroutine read_angles(kanal,*)
1167       implicit real*8 (a-h,o-z)
1168       include 'DIMENSIONS'
1169       include 'COMMON.GEO'
1170       include 'COMMON.VAR'
1171       include 'COMMON.CHAIN'
1172       include 'COMMON.IOUNITS'
1173       include 'COMMON.CONTROL'
1174 c Read angles from input
1175 c
1176        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1177        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1178        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1179        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1180
1181        do i=1,nres
1182 c 9/7/01 avoid 180 deg valence angle
1183         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1184 c
1185         theta(i)=deg2rad*theta(i)
1186         phi(i)=deg2rad*phi(i)
1187         alph(i)=deg2rad*alph(i)
1188         omeg(i)=deg2rad*omeg(i)
1189        enddo
1190       return
1191    10 return1
1192       end
1193 c----------------------------------------------------------------------------
1194       subroutine reada(rekord,lancuch,wartosc,default)
1195       implicit none
1196       character*(*) rekord,lancuch
1197       double precision wartosc,default
1198       integer ilen,iread
1199       external ilen
1200       iread=index(rekord,lancuch)
1201       if (iread.eq.0) then
1202         wartosc=default
1203         return
1204       endif
1205       iread=iread+ilen(lancuch)+1
1206       read (rekord(iread:),*,err=10,end=10) wartosc
1207       return
1208   10  wartosc=default
1209       return
1210       end
1211 c----------------------------------------------------------------------------
1212       subroutine readi(rekord,lancuch,wartosc,default)
1213       implicit none
1214       character*(*) rekord,lancuch
1215       integer wartosc,default
1216       integer ilen,iread
1217       external ilen
1218       iread=index(rekord,lancuch)
1219       if (iread.eq.0) then
1220         wartosc=default
1221         return
1222       endif
1223       iread=iread+ilen(lancuch)+1
1224       read (rekord(iread:),*,err=10,end=10) wartosc
1225       return
1226   10  wartosc=default
1227       return
1228       end
1229 c----------------------------------------------------------------------------
1230       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1231       implicit none
1232       integer dim,i
1233       integer tablica(dim),default
1234       character*(*) rekord,lancuch
1235       character*80 aux
1236       integer ilen,iread
1237       external ilen
1238       do i=1,dim
1239         tablica(i)=default
1240       enddo
1241       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1242       if (iread.eq.0) return
1243       iread=iread+ilen(lancuch)+1
1244       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1245    10 return
1246       end
1247 c----------------------------------------------------------------------------
1248       subroutine multreada(rekord,lancuch,tablica,dim,default)
1249       implicit none
1250       integer dim,i
1251       double precision tablica(dim),default
1252       character*(*) rekord,lancuch
1253       character*80 aux
1254       integer ilen,iread
1255       external ilen
1256       do i=1,dim
1257         tablica(i)=default
1258       enddo
1259       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1260       if (iread.eq.0) return
1261       iread=iread+ilen(lancuch)+1
1262       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1263    10 return
1264       end
1265 c----------------------------------------------------------------------------
1266       subroutine openunits
1267       implicit real*8 (a-h,o-z)
1268       include 'DIMENSIONS'
1269 #ifdef MPI
1270       include 'mpif.h'
1271       character*16 form,nodename
1272       integer nodelen
1273 #endif
1274       include 'COMMON.SETUP'
1275       include 'COMMON.IOUNITS'
1276 c      include 'COMMON.MD'
1277       include 'COMMON.CONTROL'
1278       integer lenpre,lenpot,ilen,lentmp
1279       external ilen
1280       character*3 out1file_text,ucase
1281       character*3 ll
1282       external ucase
1283 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1284       call getenv_loc("PREFIX",prefix)
1285       pref_orig = prefix
1286       call getenv_loc("POT",pot)
1287       call getenv_loc("DIRTMP",tmpdir)
1288       call getenv_loc("CURDIR",curdir)
1289       call getenv_loc("OUT1FILE",out1file_text)
1290 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1291       out1file_text=ucase(out1file_text)
1292       if (out1file_text(1:1).eq."Y") then
1293         out1file=.true.
1294       else
1295         out1file=fg_rank.gt.0
1296       endif
1297       lenpre=ilen(prefix)
1298       lenpot=ilen(pot)
1299       lentmp=ilen(tmpdir)
1300       if (lentmp.gt.0) then
1301           write (*,'(80(1h!))')
1302           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1303           write (*,'(80(1h!))')
1304           write (*,*)"All output files will be on node /tmp directory."
1305 #ifdef MPI
1306         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1307         if (me.eq.king) then
1308           write (*,*) "The master node is ",nodename
1309         else if (fg_rank.eq.0) then
1310           write (*,*) "I am the CG slave node ",nodename
1311         else
1312           write (*,*) "I am the FG slave node ",nodename
1313         endif
1314 #endif
1315         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1316         lenpre = lentmp+lenpre+1
1317       endif
1318       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1319 C Get the names and open the input files
1320 #if defined(WINIFL) || defined(WINPGI)
1321       open(1,file=pref_orig(:ilen(pref_orig))//
1322      &  '.inp',status='old',readonly,shared)
1323        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1324 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1325 C Get parameter filenames and open the parameter files.
1326       call getenv_loc('BONDPAR',bondname)
1327       open (ibond,file=bondname,status='old',readonly,shared)
1328       call getenv_loc('THETPAR',thetname)
1329       open (ithep,file=thetname,status='old',readonly,shared)
1330 #ifndef CRYST_THETA
1331       call getenv_loc('THETPARPDB',thetname_pdb)
1332       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1333 #endif
1334       call getenv_loc('ROTPAR',rotname)
1335       open (irotam,file=rotname,status='old',readonly,shared)
1336 #ifndef CRYST_SC
1337       call getenv_loc('ROTPARPDB',rotname_pdb)
1338       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1339 #endif
1340       call getenv_loc('TORPAR',torname)
1341       open (itorp,file=torname,status='old',readonly,shared)
1342       call getenv_loc('TORDPAR',tordname)
1343       open (itordp,file=tordname,status='old',readonly,shared)
1344       call getenv_loc('FOURIER',fouriername)
1345       open (ifourier,file=fouriername,status='old',readonly,shared)
1346       call getenv_loc('ELEPAR',elename)
1347       open (ielep,file=elename,status='old',readonly,shared)
1348       call getenv_loc('SIDEPAR',sidename)
1349       open (isidep,file=sidename,status='old',readonly,shared)
1350 #elif (defined CRAY) || (defined AIX)
1351       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1352      &  action='read')
1353 c      print *,"Processor",myrank," opened file 1"
1354       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1355 c      print *,"Processor",myrank," opened file 9"
1356 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1357 C Get parameter filenames and open the parameter files.
1358       call getenv_loc('BONDPAR',bondname)
1359       open (ibond,file=bondname,status='old',action='read')
1360 c      print *,"Processor",myrank," opened file IBOND"
1361       call getenv_loc('THETPAR',thetname)
1362       open (ithep,file=thetname,status='old',action='read')
1363 c      print *,"Processor",myrank," opened file ITHEP"
1364 #ifndef CRYST_THETA
1365       call getenv_loc('THETPARPDB',thetname_pdb)
1366       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1367 #endif
1368       call getenv_loc('ROTPAR',rotname)
1369       open (irotam,file=rotname,status='old',action='read')
1370 c      print *,"Processor",myrank," opened file IROTAM"
1371 #ifndef CRYST_SC
1372       call getenv_loc('ROTPARPDB',rotname_pdb)
1373       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1374 #endif
1375       call getenv_loc('TORPAR',torname)
1376       open (itorp,file=torname,status='old',action='read')
1377 c      print *,"Processor",myrank," opened file ITORP"
1378       call getenv_loc('TORDPAR',tordname)
1379       open (itordp,file=tordname,status='old',action='read')
1380 c      print *,"Processor",myrank," opened file ITORDP"
1381       call getenv_loc('SCCORPAR',sccorname)
1382       open (isccor,file=sccorname,status='old',action='read')
1383 c      print *,"Processor",myrank," opened file ISCCOR"
1384       call getenv_loc('FOURIER',fouriername)
1385       open (ifourier,file=fouriername,status='old',action='read')
1386 c      print *,"Processor",myrank," opened file IFOURIER"
1387       call getenv_loc('ELEPAR',elename)
1388       open (ielep,file=elename,status='old',action='read')
1389 c      print *,"Processor",myrank," opened file IELEP"
1390       call getenv_loc('SIDEPAR',sidename)
1391       open (isidep,file=sidename,status='old',action='read')
1392 c      print *,"Processor",myrank," opened file ISIDEP"
1393 c      print *,"Processor",myrank," opened parameter files"
1394 #elif (defined G77)
1395       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1396       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1397 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1398 C Get parameter filenames and open the parameter files.
1399       call getenv_loc('BONDPAR',bondname)
1400       open (ibond,file=bondname,status='old')
1401       call getenv_loc('THETPAR',thetname)
1402       open (ithep,file=thetname,status='old')
1403 #ifndef CRYST_THETA
1404       call getenv_loc('THETPARPDB',thetname_pdb)
1405       open (ithep_pdb,file=thetname_pdb,status='old')
1406 #endif
1407       call getenv_loc('ROTPAR',rotname)
1408       open (irotam,file=rotname,status='old')
1409 #ifndef CRYST_SC
1410       call getenv_loc('ROTPARPDB',rotname_pdb)
1411       open (irotam_pdb,file=rotname_pdb,status='old')
1412 #endif
1413       call getenv_loc('TORPAR',torname)
1414       open (itorp,file=torname,status='old')
1415       call getenv_loc('TORDPAR',tordname)
1416       open (itordp,file=tordname,status='old')
1417       call getenv_loc('SCCORPAR',sccorname)
1418       open (isccor,file=sccorname,status='old')
1419       call getenv_loc('FOURIER',fouriername)
1420       open (ifourier,file=fouriername,status='old')
1421       call getenv_loc('ELEPAR',elename)
1422       open (ielep,file=elename,status='old')
1423       call getenv_loc('SIDEPAR',sidename)
1424       open (isidep,file=sidename,status='old')
1425 #else
1426       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1427      &  readonly)
1428        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1429 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1430 C Get parameter filenames and open the parameter files.
1431       call getenv_loc('BONDPAR',bondname)
1432       open (ibond,file=bondname,status='old',readonly)
1433       call getenv_loc('THETPAR',thetname)
1434       open (ithep,file=thetname,status='old',readonly)
1435 #ifndef CRYST_THETA
1436       call getenv_loc('THETPARPDB',thetname_pdb)
1437       print *,"thetname_pdb ",thetname_pdb
1438       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1439       print *,ithep_pdb," opened"
1440 #endif
1441       call getenv_loc('ROTPAR',rotname)
1442       open (irotam,file=rotname,status='old',readonly)
1443 #ifndef CRYST_SC
1444       call getenv_loc('ROTPARPDB',rotname_pdb)
1445       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1446 #endif
1447       call getenv_loc('TORPAR',torname)
1448       open (itorp,file=torname,status='old',readonly)
1449       call getenv_loc('TORDPAR',tordname)
1450       open (itordp,file=tordname,status='old',readonly)
1451       call getenv_loc('SCCORPAR',sccorname)
1452       open (isccor,file=sccorname,status='old',readonly)
1453       call getenv_loc('FOURIER',fouriername)
1454       open (ifourier,file=fouriername,status='old',readonly)
1455       call getenv_loc('ELEPAR',elename)
1456       open (ielep,file=elename,status='old',readonly)
1457       call getenv_loc('SIDEPAR',sidename)
1458       open (isidep,file=sidename,status='old',readonly)
1459 #endif
1460 #ifndef OLDSCP
1461 C
1462 C 8/9/01 In the newest version SCp interaction constants are read from a file
1463 C Use -DOLDSCP to use hard-coded constants instead.
1464 C
1465       call getenv_loc('SCPPAR',scpname)
1466 #if defined(WINIFL) || defined(WINPGI)
1467       open (iscpp,file=scpname,status='old',readonly,shared)
1468 #elif (defined CRAY)  || (defined AIX)
1469       open (iscpp,file=scpname,status='old',action='read')
1470 #elif (defined G77)
1471       open (iscpp,file=scpname,status='old')
1472 #else
1473       open (iscpp,file=scpname,status='old',readonly)
1474 #endif
1475 #endif
1476       call getenv_loc('PATTERN',patname)
1477 #if defined(WINIFL) || defined(WINPGI)
1478       open (icbase,file=patname,status='old',readonly,shared)
1479 #elif (defined CRAY)  || (defined AIX)
1480       open (icbase,file=patname,status='old',action='read')
1481 #elif (defined G77)
1482       open (icbase,file=patname,status='old')
1483 #else
1484       open (icbase,file=patname,status='old',readonly)
1485 #endif
1486 #ifdef MPI
1487 C Open output file only for CG processes
1488 c      print *,"Processor",myrank," fg_rank",fg_rank
1489       if (fg_rank.eq.0) then
1490
1491       if (nodes.eq.1) then
1492         npos=3
1493       else
1494         npos = dlog10(dfloat(nodes-1))+1
1495       endif
1496       if (npos.lt.3) npos=3
1497       write (liczba,'(i1)') npos
1498       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1499      &  //')'
1500       write (liczba,form) me
1501       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1502      &  liczba(:ilen(liczba))
1503       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1504      &  //'.int'
1505       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1506      &  //'.pdb'
1507       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1508      &  liczba(:ilen(liczba))//'.mol2'
1509       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1510      &  liczba(:ilen(liczba))//'.stat'
1511       if (lentmp.gt.0)
1512      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1513      &      //liczba(:ilen(liczba))//'.stat')
1514       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1515      &  //'.rst'
1516 c      if(usampl) then
1517 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1518 c     & liczba(:ilen(liczba))//'.const'
1519 c      endif
1520
1521       endif
1522 #else
1523       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1524       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1525       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1526       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1527       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1528       if (lentmp.gt.0)
1529      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1530      &    '.stat')
1531       rest2name=prefix(:ilen(prefix))//'.rst'
1532 c      if(usampl) then
1533 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1534 c      endif
1535 #endif
1536 #if defined(AIX) || defined(PGI)
1537       if (me.eq.king .or. .not. out1file)
1538      &   open(iout,file=outname,status='unknown')
1539 #ifdef DEBUG
1540       if (fg_rank.gt.0) then
1541         write (liczba,'(i3.3)') myrank/nfgtasks
1542         write (ll,'(bz,i3.3)') fg_rank
1543         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1544      &   status='unknown')
1545       endif
1546 #endif
1547       if(me.eq.king) then
1548        open(igeom,file=intname,status='unknown',position='append')
1549        open(ipdb,file=pdbname,status='unknown')
1550        open(imol2,file=mol2name,status='unknown')
1551        open(istat,file=statname,status='unknown',position='append')
1552       else
1553 c1out       open(iout,file=outname,status='unknown')
1554       endif
1555 #else
1556       if (me.eq.king .or. .not.out1file)
1557      &    open(iout,file=outname,status='unknown')
1558 #ifdef DEBUG
1559       if (fg_rank.gt.0) then
1560         write (liczba,'(i3.3)') myrank/nfgtasks
1561         write (ll,'(bz,i3.3)') fg_rank
1562         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1563      &   status='unknown')
1564       endif
1565 #endif
1566       if(me.eq.king) then
1567        open(igeom,file=intname,status='unknown',access='append')
1568        open(ipdb,file=pdbname,status='unknown')
1569        open(imol2,file=mol2name,status='unknown')
1570        open(istat,file=statname,status='unknown',access='append')
1571       else
1572 c1out       open(iout,file=outname,status='unknown')
1573       endif
1574 #endif
1575       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1576       csa_seed=prefix(:lenpre)//'.CSA.seed'
1577       csa_history=prefix(:lenpre)//'.CSA.history'
1578       csa_bank=prefix(:lenpre)//'.CSA.bank'
1579       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1580       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1581       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1582 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1583       csa_int=prefix(:lenpre)//'.int'
1584       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1585       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1586       csa_in=prefix(:lenpre)//'.CSA.in'
1587 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1588 C Write file names
1589       if (me.eq.king)then
1590       write (iout,'(80(1h-))')
1591       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1592       write (iout,'(80(1h-))')
1593       write (iout,*) "Input file                      : ",
1594      &  pref_orig(:ilen(pref_orig))//'.inp'
1595       write (iout,*) "Output file                     : ",
1596      &  outname(:ilen(outname))
1597       write (iout,*)
1598       write (iout,*) "Sidechain potential file        : ",
1599      &  sidename(:ilen(sidename))
1600 #ifndef OLDSCP
1601       write (iout,*) "SCp potential file              : ",
1602      &  scpname(:ilen(scpname))
1603 #endif
1604       write (iout,*) "Electrostatic potential file    : ",
1605      &  elename(:ilen(elename))
1606       write (iout,*) "Cumulant coefficient file       : ",
1607      &  fouriername(:ilen(fouriername))
1608       write (iout,*) "Torsional parameter file        : ",
1609      &  torname(:ilen(torname))
1610       write (iout,*) "Double torsional parameter file : ",
1611      &  tordname(:ilen(tordname))
1612       write (iout,*) "SCCOR parameter file : ",
1613      &  sccorname(:ilen(sccorname))
1614       write (iout,*) "Bond & inertia constant file    : ",
1615      &  bondname(:ilen(bondname))
1616       write (iout,*) "Bending parameter file          : ",
1617      &  thetname(:ilen(thetname))
1618       write (iout,*) "Rotamer parameter file          : ",
1619      &  rotname(:ilen(rotname))
1620       write (iout,*) "Threading database              : ",
1621      &  patname(:ilen(patname))
1622       if (lentmp.ne.0)
1623      &write (iout,*)" DIRTMP                          : ",
1624      &  tmpdir(:lentmp)
1625       write (iout,'(80(1h-))')
1626       endif
1627       return
1628       end
1629 c----------------------------------------------------------------------------
1630       subroutine card_concat(card)
1631       implicit real*8 (a-h,o-z)
1632       include 'DIMENSIONS'
1633       include 'COMMON.IOUNITS'
1634       character*(*) card
1635       character*80 karta,ucase
1636       external ilen
1637       read (inp,'(a)') karta
1638       karta=ucase(karta)
1639       card=' '
1640       do while (karta(80:80).eq.'&')
1641         card=card(:ilen(card)+1)//karta(:79)
1642         read (inp,'(a)') karta
1643         karta=ucase(karta)
1644       enddo
1645       card=card(:ilen(card)+1)//karta
1646       return
1647       end
1648
1649       subroutine read_dist_constr
1650       implicit real*8 (a-h,o-z)
1651       include 'DIMENSIONS'
1652 #ifdef MPI
1653       include 'mpif.h'
1654 #endif
1655       include 'COMMON.SETUP'
1656       include 'COMMON.CONTROL'
1657       include 'COMMON.CHAIN'
1658       include 'COMMON.IOUNITS'
1659       include 'COMMON.SBRIDGE'
1660       integer ifrag_(2,100),ipair_(2,100)
1661       double precision wfrag_(100),wpair_(100)
1662       character*500 controlcard
1663       write (iout,*) "Calling read_dist_constr"
1664       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1665       call flush(iout)
1666       call card_concat(controlcard)
1667       call readi(controlcard,"NFRAG",nfrag_,0)
1668       call readi(controlcard,"NPAIR",npair_,0)
1669       call readi(controlcard,"NDIST",ndist_,0)
1670       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1671       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1672       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1673       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1674       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1675 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1676 c      write (iout,*) "IFRAG"
1677 c      do i=1,nfrag_
1678 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1679 c      enddo
1680 c      write (iout,*) "IPAIR"
1681 c      do i=1,npair_
1682 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1683 c      enddo
1684       call flush(iout)
1685       do i=1,nfrag_
1686         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1687         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1688      &    ifrag_(2,i)=nstart_sup+nsup-1
1689 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1690         call flush(iout)
1691         if (wfrag_(i).gt.0.0d0) then
1692         do j=ifrag_(1,i),ifrag_(2,i)-1
1693           do k=j+1,ifrag_(2,i)
1694             write (iout,*) "j",j," k",k
1695             ddjk=dist(j,k)
1696             if (constr_dist.eq.1) then
1697             nhpb=nhpb+1
1698             ihpb(nhpb)=j
1699             jhpb(nhpb)=k
1700               dhpb(nhpb)=ddjk
1701             forcon(nhpb)=wfrag_(i)
1702             else if (constr_dist.eq.2) then
1703               if (ddjk.le.dist_cut) then
1704                 nhpb=nhpb+1
1705                 ihpb(nhpb)=j
1706                 jhpb(nhpb)=k
1707                 dhpb(nhpb)=ddjk
1708                 forcon(nhpb)=wfrag_(i)
1709               endif
1710             else
1711               nhpb=nhpb+1
1712               ihpb(nhpb)=j
1713               jhpb(nhpb)=k
1714               dhpb(nhpb)=ddjk
1715               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1716             endif
1717 #ifdef MPI
1718             if (.not.out1file .or. me.eq.king)
1719      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1720      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1721 #else
1722             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1723      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1724 #endif
1725           enddo
1726         enddo
1727         endif
1728       enddo
1729       do i=1,npair_
1730         if (wpair_(i).gt.0.0d0) then
1731         ii = ipair_(1,i)
1732         jj = ipair_(2,i)
1733         if (ii.gt.jj) then
1734           itemp=ii
1735           ii=jj
1736           jj=itemp
1737         endif
1738         do j=ifrag_(1,ii),ifrag_(2,ii)
1739           do k=ifrag_(1,jj),ifrag_(2,jj)
1740             nhpb=nhpb+1
1741             ihpb(nhpb)=j
1742             jhpb(nhpb)=k
1743             forcon(nhpb)=wpair_(i)
1744             dhpb(nhpb)=dist(j,k)
1745 #ifdef MPI
1746             if (.not.out1file .or. me.eq.king)
1747      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1748      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1749 #else
1750             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1751      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1752 #endif
1753           enddo
1754         enddo
1755         endif
1756       enddo
1757       do i=1,ndist_
1758         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1759         if (forcon(nhpb+1).gt.0.0d0) then
1760           nhpb=nhpb+1
1761           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1762 #ifdef MPI
1763           if (.not.out1file .or. me.eq.king)
1764      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1765      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1766 #else
1767           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1768      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1769 #endif
1770         endif
1771       enddo
1772       call flush(iout)
1773       return
1774       end
1775 c-------------------------------------------------------------------------------
1776 #ifdef WINIFL
1777       subroutine flush(iu)
1778       return
1779       end
1780 #endif
1781 #ifdef AIX
1782       subroutine flush(iu)
1783       call flush_(iu)
1784       return
1785       end
1786 #endif
1787 c------------------------------------------------------------------------------
1788       subroutine copy_to_tmp(source)
1789       include "DIMENSIONS"
1790       include "COMMON.IOUNITS"
1791       character*(*) source
1792       character* 256 tmpfile
1793       integer ilen
1794       external ilen
1795       logical ex
1796       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1797       inquire(file=tmpfile,exist=ex)
1798       if (ex) then
1799         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1800      &   " to temporary directory..."
1801         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1802         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1803       endif
1804       return
1805       end
1806 c------------------------------------------------------------------------------
1807       subroutine move_from_tmp(source)
1808       include "DIMENSIONS"
1809       include "COMMON.IOUNITS"
1810       character*(*) source
1811       integer ilen
1812       external ilen
1813       write (*,*) "Moving ",source(:ilen(source)),
1814      & " from temporary directory to working directory"
1815       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1816       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1817       return
1818       end
1819 c------------------------------------------------------------------------------
1820       subroutine random_init(seed)
1821 C
1822 C Initialize random number generator
1823 C
1824       implicit real*8 (a-h,o-z)
1825       include 'DIMENSIONS'
1826 #ifdef AMD64
1827       integer*8 iseedi8
1828 #endif
1829 #ifdef MPI
1830       include 'mpif.h'
1831       logical OKRandom, prng_restart
1832       real*8  r1
1833       integer iseed_array(4)
1834 #endif
1835       include 'COMMON.IOUNITS'
1836       include 'COMMON.TIME1'
1837 c      include 'COMMON.THREAD'
1838       include 'COMMON.SBRIDGE'
1839       include 'COMMON.CONTROL'
1840       include 'COMMON.MCM'
1841 c      include 'COMMON.MAP'
1842       include 'COMMON.HEADER'
1843 c      include 'COMMON.CSA'
1844       include 'COMMON.CHAIN'
1845 c      include 'COMMON.MUCA'
1846 c      include 'COMMON.MD'
1847       include 'COMMON.FFIELD'
1848       include 'COMMON.SETUP'
1849       iseed=-dint(dabs(seed))
1850       if (iseed.eq.0) then
1851         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1852      &    'Random seed undefined. The program will stop.'
1853         write (*,'(/80(1h*)/20x,a/80(1h*))')
1854      &    'Random seed undefined. The program will stop.'
1855 #ifdef MPI
1856         call mpi_finalize(mpi_comm_world,ierr)
1857 #endif
1858         stop 'Bad random seed.'
1859       endif
1860 #ifdef MPI
1861       if (fg_rank.eq.0) then
1862       seed=seed*(me+1)+1
1863 #ifdef AMD64
1864       iseedi8=dint(seed)
1865       if(me.eq.king .or. .not. out1file)
1866      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1867       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1868       OKRandom = prng_restart(me,iseedi8)
1869 #else
1870       do i=1,4
1871        tmp=65536.0d0**(4-i)
1872        iseed_array(i) = dint(seed/tmp)
1873        seed=seed-iseed_array(i)*tmp
1874       enddo
1875       if(me.eq.king .or. .not. out1file)
1876      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1877      &                 (iseed_array(i),i=1,4)
1878       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1879      &                 (iseed_array(i),i=1,4)
1880       OKRandom = prng_restart(me,iseed_array)
1881 #endif
1882       if (OKRandom) then
1883         r1=ran_number(0.0D0,1.0D0)
1884         if(me.eq.king .or. .not. out1file)
1885      &   write (iout,*) 'ran_num',r1
1886         if (r1.lt.0.0d0) OKRandom=.false.
1887       endif
1888       if (.not.OKRandom) then
1889         write (iout,*) 'PRNG IS NOT WORKING!!!'
1890         print *,'PRNG IS NOT WORKING!!!'
1891         if (me.eq.0) then
1892          call flush(iout)
1893          call mpi_abort(mpi_comm_world,error_msg,ierr)
1894          stop
1895         else
1896          write (iout,*) 'too many processors for parallel prng'
1897          write (*,*) 'too many processors for parallel prng'
1898          call flush(iout)
1899          stop
1900         endif
1901       endif
1902       endif
1903 #else
1904       call vrndst(iseed)
1905       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1906 #endif
1907       return
1908       end