2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
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 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 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
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
52 c print *,"Processor",myrank," leaves READRTNS"
55 C-------------------------------------------------------------------------------
56 subroutine read_control
60 implicit real*8 (a-h,o-z)
64 logical OKRandom, prng_restart
67 include 'COMMON.IOUNITS'
68 include 'COMMON.TIME1'
69 include 'COMMON.THREAD'
70 include 'COMMON.SBRIDGE'
71 include 'COMMON.CONTROL'
74 include 'COMMON.HEADER'
76 include 'COMMON.CHAIN'
79 include 'COMMON.FFIELD'
80 include 'COMMON.INTERACT'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
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 readi(controlcard,'SYM',symetr,1)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
181 crc else if (index(controlcard,'ZSCORE').gt.0) then
183 crc ZSCORE is rm from UNRES, modecalc=9 is available
186 cfcm else if (index(controlcard,'MCMF').gt.0) then
188 else if (index(controlcard,'SOFTREG').gt.0) then
190 else if (index(controlcard,'CHECK_BOND').gt.0) then
192 else if (index(controlcard,'TEST').gt.0) then
194 else if (index(controlcard,'MD').gt.0) then
196 else if (index(controlcard,'RE ').gt.0) then
200 lmuca=index(controlcard,'MUCA').gt.0
201 call readi(controlcard,'MUCADYN',mucadyn,0)
202 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
203 if (lmuca .and. (me.eq.king .or. .not.out1file ))
205 write (iout,*) 'MUCADYN=',mucadyn
206 write (iout,*) 'MUCASMOOTH=',muca_smooth
209 iscode=index(controlcard,'ONE_LETTER')
210 indphi=index(controlcard,'PHI')
211 indback=index(controlcard,'BACK')
212 iranconf=index(controlcard,'RAND_CONF')
213 i2ndstr=index(controlcard,'USE_SEC_PRED')
214 gradout=index(controlcard,'GRADOUT').gt.0
215 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
217 if (me.eq.king .or. .not.out1file )
218 & write (iout,*) "DISTCHAINMAX",distchainmax
220 if(me.eq.king.or..not.out1file)
221 & write (iout,'(2a)') diagmeth(kdiag),
222 & ' routine used to diagonalize matrices.'
225 c--------------------------------------------------------------------------
226 subroutine read_REMDpar
230 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
233 include 'COMMON.TIME1'
236 include 'COMMON.LANGEVIN'
238 include 'COMMON.LANGEVIN.lang0'
240 include 'COMMON.INTERACT'
241 include 'COMMON.NAMES'
243 include 'COMMON.REMD'
244 include 'COMMON.CONTROL'
245 include 'COMMON.SETUP'
247 character*320 controlcard
248 character*3200 controlcard1
249 integer iremd_m_total
251 if(me.eq.king.or..not.out1file)
252 & write (iout,*) "REMD setup"
254 call card_concat(controlcard)
255 call readi(controlcard,"NREP",nrep,3)
256 call readi(controlcard,"NSTEX",nstex,1000)
257 call reada(controlcard,"RETMIN",retmin,10.0d0)
258 call reada(controlcard,"RETMAX",retmax,1000.0d0)
259 mremdsync=(index(controlcard,'SYNC').gt.0)
260 call readi(controlcard,"NSYN",i_sync_step,100)
261 restart1file=(index(controlcard,'REST1FILE').gt.0)
262 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
263 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
264 if(max_cache_traj_use.gt.max_cache_traj)
265 & max_cache_traj_use=max_cache_traj
266 if(me.eq.king.or..not.out1file) then
267 cd if (traj1file) then
268 crc caching is in testing - NTWX is not ignored
269 cd write (iout,*) "NTWX value is ignored"
270 cd write (iout,*) " trajectory is stored to one file by master"
271 cd write (iout,*) " before exchange at NSTEX intervals"
273 write (iout,*) "NREP= ",nrep
274 write (iout,*) "NSTEX= ",nstex
275 write (iout,*) "SYNC= ",mremdsync
276 write (iout,*) "NSYN= ",i_sync_step
277 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
280 if (index(controlcard,'TLIST').gt.0) then
282 call card_concat(controlcard1)
283 read(controlcard1,*) (remd_t(i),i=1,nrep)
284 if(me.eq.king.or..not.out1file)
285 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
288 if (index(controlcard,'MLIST').gt.0) then
290 call card_concat(controlcard1)
291 read(controlcard1,*) (remd_m(i),i=1,nrep)
292 if(me.eq.king.or..not.out1file) then
293 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
296 iremd_m_total=iremd_m_total+remd_m(i)
298 write (iout,*) 'Total number of replicas ',iremd_m_total
301 if(me.eq.king.or..not.out1file)
302 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
305 c--------------------------------------------------------------------------
306 subroutine read_MDpar
310 implicit real*8 (a-h,o-z)
312 include 'COMMON.IOUNITS'
313 include 'COMMON.TIME1'
316 include 'COMMON.LANGEVIN'
318 include 'COMMON.LANGEVIN.lang0'
320 include 'COMMON.INTERACT'
321 include 'COMMON.NAMES'
323 include 'COMMON.SETUP'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SPLITELE'
327 character*320 controlcard
329 call card_concat(controlcard)
330 call readi(controlcard,"NSTEP",n_timestep,1000000)
331 call readi(controlcard,"NTWE",ntwe,100)
332 call readi(controlcard,"NTWX",ntwx,1000)
333 call reada(controlcard,"DT",d_time,1.0d-1)
334 call reada(controlcard,"DVMAX",dvmax,2.0d1)
335 call reada(controlcard,"DAMAX",damax,1.0d1)
336 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
337 call readi(controlcard,"LANG",lang,0)
338 RESPA = index(controlcard,"RESPA") .gt. 0
339 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
340 ntime_split0=ntime_split
341 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
342 ntime_split0=ntime_split
343 call reada(controlcard,"R_CUT",r_cut,2.0d0)
344 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
345 rest = index(controlcard,"REST").gt.0
346 tbf = index(controlcard,"TBF").gt.0
347 usampl = index(controlcard,"USAMPL").gt.0
349 mdpdb = index(controlcard,"MDPDB").gt.0
350 call reada(controlcard,"T_BATH",t_bath,300.0d0)
351 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
352 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
353 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
354 if (count_reset_moment.eq.0) count_reset_moment=1000000000
355 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
356 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
357 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
358 if (count_reset_vel.eq.0) count_reset_vel=1000000000
359 large = index(controlcard,"LARGE").gt.0
360 print_compon = index(controlcard,"PRINT_COMPON").gt.0
361 rattle = index(controlcard,"RATTLE").gt.0
362 c if performing umbrella sampling, fragments constrained are read from the fragment file
368 if(me.eq.king.or..not.out1file) then
370 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
372 write (iout,'(a)') "The units are:"
373 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
374 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
375 & " acceleration: angstrom/(48.9 fs)**2"
376 write (iout,'(a)') "energy: kcal/mol, temperature: K"
378 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
379 write (iout,'(a60,f10.5,a)')
380 & "Initial time step of numerical integration:",d_time,
382 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
384 write (iout,'(2a,i4,a)')
385 & "A-MTS algorithm used; initial time step for fast-varying",
386 & " short-range forces split into",ntime_split," steps."
387 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
388 & r_cut," lambda",rlamb
390 write (iout,'(2a,f10.5)')
391 & "Maximum acceleration threshold to reduce the time step",
392 & "/increase split number:",damax
393 write (iout,'(2a,f10.5)')
394 & "Maximum predicted energy drift to reduce the timestep",
395 & "/increase split number:",edriftmax
396 write (iout,'(a60,f10.5)')
397 & "Maximum velocity threshold to reduce velocities:",dvmax
398 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
399 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
400 if (rattle) write (iout,'(a60)')
401 & "Rattle algorithm used to constrain the virtual bonds"
405 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
406 call reada(controlcard,"RWAT",rwat,1.4d0)
407 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
408 surfarea=index(controlcard,"SURFAREA").gt.0
409 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
410 if(me.eq.king.or..not.out1file)then
411 write (iout,'(/a,$)') "Langevin dynamics calculation"
414 & " with direct integration of Langevin equations"
415 else if (lang.eq.2) then
416 write (iout,'(a/)') " with TINKER stochasic MD integrator"
417 else if (lang.eq.3) then
418 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
419 else if (lang.eq.4) then
420 write (iout,'(a/)') " in overdamped mode"
422 write (iout,'(//a,i5)')
423 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
426 write (iout,'(a60,f10.5)') "Temperature:",t_bath
427 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
428 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
429 write (iout,'(a60,f10.5)')
430 & "Scaling factor of the friction forces:",scal_fric
431 if (surfarea) write (iout,'(2a,i10,a)')
432 & "Friction coefficients will be scaled by solvent-accessible",
433 & " surface area every",reset_fricmat," steps."
435 c Calculate friction coefficients and bounds of stochastic forces
436 eta=6*pi*cPoise*etawat
437 if(me.eq.king.or..not.out1file)
438 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
440 gamp=scal_fric*(pstok+rwat)*eta
441 stdfp=dsqrt(2*Rb*t_bath/d_time)
443 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
444 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
446 if(me.eq.king.or..not.out1file)then
447 write (iout,'(/2a/)')
448 & "Radii of site types and friction coefficients and std's of",
449 & " stochastic forces of fully exposed sites"
450 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
452 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
453 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
457 if(me.eq.king.or..not.out1file)then
458 write (iout,'(a)') "Berendsen bath calculation"
459 write (iout,'(a60,f10.5)') "Temperature:",t_bath
460 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
462 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
463 & count_reset_moment," steps"
465 & write (iout,'(a,i10,a)')
466 & "Velocities will be reset at random every",count_reset_vel,
470 if(me.eq.king.or..not.out1file)
471 & write (iout,'(a31)') "Microcanonical mode calculation"
473 if(me.eq.king.or..not.out1file)then
474 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
476 write(iout,*) "MD running with constraints."
477 write(iout,*) "Equilibration time ", eq_time, " mtus."
478 write(iout,*) "Constraining ", nfrag," fragments."
479 write(iout,*) "Length of each fragment, weight and q0:"
481 write (iout,*) "Set of restraints #",iset
483 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
484 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
486 write(iout,*) "constraints between ", npair, "fragments."
487 write(iout,*) "constraint pairs, weights and q0:"
489 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
490 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
492 write(iout,*) "angle constraints within ", nfrag_back,
493 & "backbone fragments."
494 write(iout,*) "fragment, weights:"
496 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
497 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
498 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
501 iset=mod(kolor,nset)+1
504 if(me.eq.king.or..not.out1file)
505 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
508 c------------------------------------------------------------------------------
511 C Read molecular data.
513 implicit real*8 (a-h,o-z)
519 include 'COMMON.IOUNITS'
522 include 'COMMON.INTERACT'
523 include 'COMMON.LOCAL'
524 include 'COMMON.NAMES'
525 include 'COMMON.CHAIN'
526 include 'COMMON.FFIELD'
527 include 'COMMON.SBRIDGE'
528 include 'COMMON.HEADER'
529 include 'COMMON.CONTROL'
530 include 'COMMON.DBASE'
531 include 'COMMON.THREAD'
532 include 'COMMON.CONTACTS'
533 include 'COMMON.TORCNSTR'
534 include 'COMMON.TIME1'
535 include 'COMMON.BOUNDS'
537 include 'COMMON.SETUP'
538 character*4 sequence(maxres)
540 double precision x(maxvar)
541 character*256 pdbfile
542 character*400 weightcard
543 character*80 weightcard_t,ucase
544 dimension itype_pdb(maxres)
545 common /pizda/ itype_pdb
546 logical seq_comp,fail
547 double precision energia(0:n_ene)
553 C Read weights of the subsequent energy terms.
554 call card_concat(weightcard)
555 call reada(weightcard,'WLONG',wlong,1.0D0)
556 call reada(weightcard,'WSC',wsc,wlong)
557 call reada(weightcard,'WSCP',wscp,wlong)
558 call reada(weightcard,'WELEC',welec,1.0D0)
559 call reada(weightcard,'WVDWPP',wvdwpp,welec)
560 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
561 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
562 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
563 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
564 call reada(weightcard,'WTURN3',wturn3,1.0D0)
565 call reada(weightcard,'WTURN4',wturn4,1.0D0)
566 call reada(weightcard,'WTURN6',wturn6,1.0D0)
567 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
568 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
569 call reada(weightcard,'WBOND',wbond,1.0D0)
570 call reada(weightcard,'WTOR',wtor,1.0D0)
571 call reada(weightcard,'WTORD',wtor_d,1.0D0)
572 call reada(weightcard,'WANG',wang,1.0D0)
573 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
574 call reada(weightcard,'SCAL14',scal14,0.4D0)
575 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
576 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
577 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
578 call reada(weightcard,'TEMP0',temp0,300.0d0)
579 if (index(weightcard,'SOFT').gt.0) ipot=6
580 C 12/1/95 Added weight for the multi-body term WCORR
581 call reada(weightcard,'WCORRH',wcorr,1.0D0)
582 if (wcorr4.gt.0.0d0) wcorr=wcorr4
602 if(me.eq.king.or..not.out1file)
603 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
604 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
606 10 format (/'Energy-term weights (unscaled):'//
607 & 'WSCC= ',f10.6,' (SC-SC)'/
608 & 'WSCP= ',f10.6,' (SC-p)'/
609 & 'WELEC= ',f10.6,' (p-p electr)'/
610 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
611 & 'WBOND= ',f10.6,' (stretching)'/
612 & 'WANG= ',f10.6,' (bending)'/
613 & 'WSCLOC= ',f10.6,' (SC local)'/
614 & 'WTOR= ',f10.6,' (torsional)'/
615 & 'WTORD= ',f10.6,' (double torsional)'/
616 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
617 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
618 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
619 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
620 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
621 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
622 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
623 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
624 & 'WTURN6= ',f10.6,' (turns, 6th order)')
625 if(me.eq.king.or..not.out1file)then
626 if (wcorr4.gt.0.0d0) then
627 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
628 & 'between contact pairs of peptide groups'
629 write (iout,'(2(a,f5.3/))')
630 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
631 & 'Range of quenching the correlation terms:',2*delt_corr
632 else if (wcorr.gt.0.0d0) then
633 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
634 & 'between contact pairs of peptide groups'
636 write (iout,'(a,f8.3)')
637 & 'Scaling factor of 1,4 SC-p interactions:',scal14
638 write (iout,'(a,f8.3)')
639 & 'General scaling factor of SC-p interactions:',scalscp
641 r0_corr=cutoff_corr-delt_corr
643 aad(i,1)=scalscp*aad(i,1)
644 aad(i,2)=scalscp*aad(i,2)
645 bad(i,1)=scalscp*bad(i,1)
646 bad(i,2)=scalscp*bad(i,2)
648 call rescale_weights(t_bath)
649 if(me.eq.king.or..not.out1file)
650 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
651 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
653 22 format (/'Energy-term weights (scaled):'//
654 & 'WSCC= ',f10.6,' (SC-SC)'/
655 & 'WSCP= ',f10.6,' (SC-p)'/
656 & 'WELEC= ',f10.6,' (p-p electr)'/
657 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
658 & 'WBOND= ',f10.6,' (stretching)'/
659 & 'WANG= ',f10.6,' (bending)'/
660 & 'WSCLOC= ',f10.6,' (SC local)'/
661 & 'WTOR= ',f10.6,' (torsional)'/
662 & 'WTORD= ',f10.6,' (double torsional)'/
663 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
664 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
665 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
666 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
667 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
668 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
669 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
670 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
671 & 'WTURN6= ',f10.6,' (turns, 6th order)')
672 if(me.eq.king.or..not.out1file)
673 & write (iout,*) "Reference temperature for weights calculation:",
675 call reada(weightcard,"D0CM",d0cm,3.78d0)
676 call reada(weightcard,"AKCM",akcm,15.1d0)
677 call reada(weightcard,"AKTH",akth,11.0d0)
678 call reada(weightcard,"AKCT",akct,12.0d0)
679 call reada(weightcard,"V1SS",v1ss,-1.08d0)
680 call reada(weightcard,"V2SS",v2ss,7.61d0)
681 call reada(weightcard,"V3SS",v3ss,13.7d0)
682 call reada(weightcard,"EBR",ebr,-5.50D0)
683 call reada(weightcard,"ATRISS",atriss,0.301D0)
684 call reada(weightcard,"BTRISS",btriss,0.021D0)
685 call reada(weightcard,"CTRISS",ctriss,1.001D0)
686 call reada(weightcard,"DTRISS",dtriss,1.001D0)
687 write (iout,*) "ATRISS=", atriss
688 write (iout,*) "BTRISS=", btriss
689 write (iout,*) "CTRISS=", ctriss
690 write (iout,*) "DTRISS=", dtriss
691 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
693 dyn_ss_mask(i)=.false.
697 dyn_ssbond_ij(i,j)=1.0d300
700 call reada(weightcard,"HT",Ht,0.0D0)
702 ss_depth=ebr/wsc-0.25*eps(1,1)
703 Ht=Ht/wsc-0.25*eps(1,1)
704 akcm=akcm*wstrain/wsc
705 akth=akth*wstrain/wsc
706 akct=akct*wstrain/wsc
707 v1ss=v1ss*wstrain/wsc
708 v2ss=v2ss*wstrain/wsc
709 v3ss=v3ss*wstrain/wsc
711 if (wstrain.ne.0.0) then
712 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
718 if(me.eq.king.or..not.out1file) then
719 write (iout,*) "Parameters of the SS-bond potential:"
720 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
722 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
723 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
724 write (iout,*)" HT",Ht
725 print *,'indpdb=',indpdb,' pdbref=',pdbref
727 if (indpdb.gt.0 .or. pdbref) then
728 read(inp,'(a)') pdbfile
729 if(me.eq.king.or..not.out1file)
730 & write (iout,'(2a)') 'PDB data will be read from file ',
731 & pdbfile(:ilen(pdbfile))
732 open(ipdbin,file=pdbfile,status='old',err=33)
734 33 write (iout,'(a)') 'Error opening PDB file.'
737 c write (iout,*) 'Begin reading pdb data'
740 c write (iout,*) 'Finished reading pdb data'
742 if(me.eq.king.or..not.out1file)
743 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
744 & ' nstart_sup=',nstart_sup
746 itype_pdb(i)=itype(i)
750 nct=nstart_sup+nsup-1
751 call contact(.false.,ncont_ref,icont_ref,co)
754 if(me.eq.king.or..not.out1file)
755 & write(iout,*)'Adding sidechains'
759 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
762 do while (fail.and.nsi.le.maxsi)
763 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
766 if(fail) write(iout,*)'Adding sidechain failed for res ',
767 & i,' after ',nsi,' trials'
772 if (indpdb.eq.0) then
773 C Read sequence if not taken from the pdb file.
775 c print *,'nres=',nres
776 if (iscode.gt.0) then
777 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
779 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
781 C Convert sequence to numeric code
783 itype(i)=rescode(i,sequence(i),iscode)
785 C Assign initial virtual bond lengths
791 vbld(i+nres)=dsc(iabs(itype(i)))
792 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
793 c write (iout,*) "i",i," itype",itype(i),
794 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
798 c print '(20i4)',(itype(i),i=1,nres)
801 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
803 if (itype(i).eq.ntyp1) then
807 else if (iabs(itype(i+1)).ne.20) then
809 else if (iabs(itype(i)).ne.20) then
816 if(me.eq.king.or..not.out1file)then
817 write (iout,*) "ITEL"
819 write (iout,*) i,itype(i),itel(i)
821 print *,'Call Read_Bridge.'
824 C 8/13/98 Set limits to generating the dihedral angles
829 read (inp,*) ndih_constr
830 if (ndih_constr.gt.0) then
832 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
833 if(me.eq.king.or..not.out1file)then
835 & 'There are',ndih_constr,' constraints on phi angles.'
837 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
841 phi0(i)=deg2rad*phi0(i)
842 drange(i)=deg2rad*drange(i)
844 if(me.eq.king.or..not.out1file)
845 & write (iout,*) 'FTORS',ftors
848 phibound(1,ii) = phi0(i)-drange(i)
849 phibound(2,ii) = phi0(i)+drange(i)
856 write (iout,'(a)') 'Boundaries in phi angle sampling:'
858 write (iout,'(a3,i5,2f10.1)')
859 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
865 cd print *,'NNT=',NNT,' NCT=',NCT
866 if (itype(1).eq.ntyp1) nnt=2
867 if (itype(nres).eq.ntyp1) nct=nct-1
869 if(me.eq.king.or..not.out1file)
870 & write (iout,'(a,i3)') 'nsup=',nsup
872 if (nsup.le.(nct-nnt+1)) then
873 do i=0,nct-nnt+1-nsup
874 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
880 & 'Error - sequences to be superposed do not match.'
883 do i=0,nsup-(nct-nnt+1)
884 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
886 nstart_sup=nstart_sup+i
892 & 'Error - sequences to be superposed do not match.'
895 if (nsup.eq.0) nsup=nct-nnt
896 if (nstart_sup.eq.0) nstart_sup=nnt
897 if (nstart_seq.eq.0) nstart_seq=nnt
898 if(me.eq.king.or..not.out1file)
899 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
900 & ' nstart_seq=',nstart_seq
902 c--- Zscore rms -------
903 if (nz_start.eq.0) nz_start=nnt
904 if (nz_end.eq.0 .and. nsup.gt.0) then
906 else if (nz_end.eq.0) then
909 if(me.eq.king.or..not.out1file)then
910 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
911 write (iout,*) 'IZ_SC=',iz_sc
913 c----------------------
916 if (.not.pdbref) then
917 call read_angles(inp,*38)
919 38 write (iout,'(a)') 'Error reading reference structure.'
921 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
922 stop 'Error reading reference structure'
926 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
936 call contact(.true.,ncont_ref,icont_ref,co)
940 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
942 if (constr_dist.gt.0) call read_dist_constr
943 write (iout,*) "After read_dist_constr nhpb",nhpb
945 if(me.eq.king.or..not.out1file)
946 & write (iout,*) 'Contact order:',co
948 if(me.eq.king.or..not.out1file)
949 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
952 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
954 if(me.eq.king.or..not.out1file)
955 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
956 & icont_ref(1,i),' ',
957 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
961 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
962 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
963 & modecalc.ne.10) then
964 C If input structure hasn't been supplied from the PDB file read or generate
966 if (iranconf.eq.0 .and. .not. extconf) then
967 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
968 & write (iout,'(a)') 'Initial geometry will be read in.'
970 read(inp,'(8f10.5)',end=36,err=36)
971 & ((c(l,k),l=1,3),k=1,nres),
972 & ((c(l,k+nres),l=1,3),k=nnt,nct)
973 write (iout,*) "Exit READ_CART"
974 write (iout,'(8f10.5)')
975 & ((c(l,k),l=1,3),k=1,nres),
976 & ((c(l,k+nres),l=1,3),k=nnt,nct)
977 call int_from_cart1(.true.)
978 write (iout,*) "Finish INT_TO_CART"
981 dc(j,i)=c(j,i+1)-c(j,i)
982 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
986 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
988 dc(j,i+nres)=c(j,i+nres)-c(j,i)
989 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
995 call read_angles(inp,*36)
998 36 write (iout,'(a)') 'Error reading angle file.'
1000 call mpi_finalize( MPI_COMM_WORLD,IERR )
1002 stop 'Error reading angle file.'
1004 else if (extconf) then
1005 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1006 & write (iout,'(a)') 'Extended chain initial geometry.'
1008 theta(i)=90d0*deg2rad
1011 phi(i)=180d0*deg2rad
1014 alph(i)=110d0*deg2rad
1017 omeg(i)=-120d0*deg2rad
1018 if (itype(i).le.0) omeg(i)=-omeg(i)
1021 if(me.eq.king.or..not.out1file)
1022 & write (iout,'(a)') 'Random-generated initial geometry.'
1026 if (me.eq.king .or. fg_rank.eq.0 .and. (
1027 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1031 call gen_rand_conf(itmp,*30)
1033 30 write (iout,*) 'Failed to generate random conformation',
1034 & ', itrial=',itrial
1035 write (*,*) 'Processor:',me,
1036 & ' Failed to generate random conformation',
1046 write (iout,'(a,i3,a)') 'Processor:',me,
1047 & ' error in generating random conformation.'
1048 write (*,'(a,i3,a)') 'Processor:',me,
1049 & ' error in generating random conformation.'
1052 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1058 call gen_rand_conf(itmp,*30)
1060 30 write (iout,*) 'Failed to generate random conformation',
1061 & ', itrial=',itrial
1062 write (*,*) 'Failed to generate random conformation',
1063 & ', itrial=',itrial
1065 write (iout,'(a,i3,a)') 'Processor:',me,
1066 & ' error in generating random conformation.'
1067 write (*,'(a,i3,a)') 'Processor:',me,
1068 & ' error in generating random conformation.'
1073 elseif (modecalc.eq.4) then
1074 read (inp,'(a)') intinname
1075 open (intin,file=intinname,status='old',err=333)
1076 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1077 & write (iout,'(a)') 'intinname',intinname
1078 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1080 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1082 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1084 stop 'Error opening angle file.'
1088 C Generate distance constraints, if the PDB structure is to be regularized.
1089 if (nthread.gt.0) then
1090 call read_threadbase
1093 if (me.eq.king .or. .not. out1file)
1095 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1096 write (iout,'(/a,i3,a)')
1097 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1098 write (iout,'(20i4)') (iss(i),i=1,ns)
1100 write(iout,*)"Running with dynamic disulfide-bond formation"
1102 write (iout,'(/a/)') 'Pre-formed links are:'
1108 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1109 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1115 if (ns.gt.0.and.dyn_ss) then
1119 forcon(i-nss)=forcon(i)
1126 dyn_ss_mask(iss(i))=.true.
1129 if (i2ndstr.gt.0) call secstrp2dihc
1130 c call geom_to_var(nvar,x)
1131 c call etotal(energia(0))
1132 c call enerprint(energia(0))
1133 c call briefout(0,etot)
1135 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1136 cd write (iout,'(a)') 'Variable list:'
1137 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1139 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1140 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1141 & 'Processor',myrank,': end reading molecular data.'
1146 c--------------------------------------------------------------------------
1147 logical function seq_comp(itypea,itypeb,length)
1149 integer length,itypea(length),itypeb(length)
1152 if (itypea(i).ne.itypeb(i)) then
1160 c-----------------------------------------------------------------------------
1161 subroutine read_bridge
1162 C Read information about disulfide bridges.
1163 implicit real*8 (a-h,o-z)
1164 include 'DIMENSIONS'
1168 include 'COMMON.IOUNITS'
1169 include 'COMMON.GEO'
1170 include 'COMMON.VAR'
1171 include 'COMMON.INTERACT'
1172 include 'COMMON.LOCAL'
1173 include 'COMMON.NAMES'
1174 include 'COMMON.CHAIN'
1175 include 'COMMON.FFIELD'
1176 include 'COMMON.SBRIDGE'
1177 include 'COMMON.HEADER'
1178 include 'COMMON.CONTROL'
1179 include 'COMMON.DBASE'
1180 include 'COMMON.THREAD'
1181 include 'COMMON.TIME1'
1182 include 'COMMON.SETUP'
1183 C Read bridging residues.
1184 read (inp,*) ns,(iss(i),i=1,ns)
1186 if(me.eq.king.or..not.out1file)
1187 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1188 C Check whether the specified bridging residues are cystines.
1190 if (itype(iss(i)).ne.1) then
1191 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1192 & 'Do you REALLY think that the residue ',
1193 & restyp(itype(iss(i))),i,
1194 & ' can form a disulfide bridge?!!!'
1195 write (*,'(2a,i3,a)')
1196 & 'Do you REALLY think that the residue ',
1197 & restyp(itype(iss(i))),i,
1198 & ' can form a disulfide bridge?!!!'
1200 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1205 C Read preformed bridges.
1207 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1209 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1212 C Check if the residues involved in bridges are in the specified list of
1213 C bridging residues.
1216 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1217 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1218 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1219 & ' contains residues present in other pairs.'
1220 write (*,'(a,i3,a)') 'Disulfide pair',i,
1221 & ' contains residues present in other pairs.'
1223 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1229 if (ihpb(i).eq.iss(j)) goto 10
1231 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1234 if (jhpb(i).eq.iss(j)) goto 20
1236 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1242 ihpb(i)=ihpb(i)+nres
1243 jhpb(i)=jhpb(i)+nres
1249 c----------------------------------------------------------------------------
1250 subroutine read_x(kanal,*)
1251 implicit real*8 (a-h,o-z)
1252 include 'DIMENSIONS'
1253 include 'COMMON.GEO'
1254 include 'COMMON.VAR'
1255 include 'COMMON.CHAIN'
1256 include 'COMMON.IOUNITS'
1257 include 'COMMON.CONTROL'
1258 include 'COMMON.LOCAL'
1259 include 'COMMON.INTERACT'
1260 c Read coordinates from input
1262 read(kanal,'(8f10.5)',end=10,err=10)
1263 & ((c(l,k),l=1,3),k=1,nres),
1264 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1267 c(j,2*nres)=c(j,nres)
1269 call int_from_cart1(.false.)
1272 dc(j,i)=c(j,i+1)-c(j,i)
1273 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1277 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1279 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1280 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1288 c----------------------------------------------------------------------------
1289 subroutine read_threadbase
1290 implicit real*8 (a-h,o-z)
1291 include 'DIMENSIONS'
1292 include 'COMMON.IOUNITS'
1293 include 'COMMON.GEO'
1294 include 'COMMON.VAR'
1295 include 'COMMON.INTERACT'
1296 include 'COMMON.LOCAL'
1297 include 'COMMON.NAMES'
1298 include 'COMMON.CHAIN'
1299 include 'COMMON.FFIELD'
1300 include 'COMMON.SBRIDGE'
1301 include 'COMMON.HEADER'
1302 include 'COMMON.CONTROL'
1303 include 'COMMON.DBASE'
1304 include 'COMMON.THREAD'
1305 include 'COMMON.TIME1'
1306 C Read pattern database for threading.
1307 read (icbase,*) nseq
1309 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1310 & nres_base(2,i),nres_base(3,i)
1311 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1313 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1314 c & nres_base(2,i),nres_base(3,i)
1315 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1319 if (weidis.eq.0.0D0) weidis=0.1D0
1328 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1329 write (iout,'(a,i5)') 'nexcl: ',nexcl
1330 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1333 c------------------------------------------------------------------------------
1334 subroutine setup_var
1335 implicit real*8 (a-h,o-z)
1336 include 'DIMENSIONS'
1337 include 'COMMON.IOUNITS'
1338 include 'COMMON.GEO'
1339 include 'COMMON.VAR'
1340 include 'COMMON.INTERACT'
1341 include 'COMMON.LOCAL'
1342 include 'COMMON.NAMES'
1343 include 'COMMON.CHAIN'
1344 include 'COMMON.FFIELD'
1345 include 'COMMON.SBRIDGE'
1346 include 'COMMON.HEADER'
1347 include 'COMMON.CONTROL'
1348 include 'COMMON.DBASE'
1349 include 'COMMON.THREAD'
1350 include 'COMMON.TIME1'
1351 C Set up variable list.
1357 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1359 ialph(i,1)=nvar+nside
1363 if (indphi.gt.0) then
1365 else if (indback.gt.0) then
1370 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1373 c----------------------------------------------------------------------------
1374 subroutine gen_dist_constr
1375 C Generate CA distance constraints.
1376 implicit real*8 (a-h,o-z)
1377 include 'DIMENSIONS'
1378 include 'COMMON.IOUNITS'
1379 include 'COMMON.GEO'
1380 include 'COMMON.VAR'
1381 include 'COMMON.INTERACT'
1382 include 'COMMON.LOCAL'
1383 include 'COMMON.NAMES'
1384 include 'COMMON.CHAIN'
1385 include 'COMMON.FFIELD'
1386 include 'COMMON.SBRIDGE'
1387 include 'COMMON.HEADER'
1388 include 'COMMON.CONTROL'
1389 include 'COMMON.DBASE'
1390 include 'COMMON.THREAD'
1391 include 'COMMON.TIME1'
1392 dimension itype_pdb(maxres)
1393 common /pizda/ itype_pdb
1395 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1396 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1397 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1399 do i=nstart_sup,nstart_sup+nsup-1
1400 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1401 cd & ' seq_pdb', restyp(itype_pdb(i))
1402 do j=i+2,nstart_sup+nsup-1
1404 ihpb(nhpb)=i+nstart_seq-nstart_sup
1405 jhpb(nhpb)=j+nstart_seq-nstart_sup
1407 dhpb(nhpb)=dist(i,j)
1410 cd write (iout,'(a)') 'Distance constraints:'
1415 cd if (ii.gt.nres) then
1420 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1421 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1422 cd & dhpb(i),forcon(i)
1426 c----------------------------------------------------------------------------
1428 implicit real*8 (a-h,o-z)
1429 include 'DIMENSIONS'
1430 include 'COMMON.MAP'
1431 include 'COMMON.IOUNITS'
1432 character*3 angid(4) /'THE','PHI','ALP','OME'/
1433 character*80 mapcard,ucase
1435 read (inp,'(a)') mapcard
1436 mapcard=ucase(mapcard)
1437 if (index(mapcard,'PHI').gt.0) then
1439 else if (index(mapcard,'THE').gt.0) then
1441 else if (index(mapcard,'ALP').gt.0) then
1443 else if (index(mapcard,'OME').gt.0) then
1446 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1447 stop 'Error - illegal variable spec in MAP card.'
1449 call readi (mapcard,'RES1',res1(imap),0)
1450 call readi (mapcard,'RES2',res2(imap),0)
1451 if (res1(imap).eq.0) then
1452 res1(imap)=res2(imap)
1453 else if (res2(imap).eq.0) then
1454 res2(imap)=res1(imap)
1456 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1458 & 'Error - illegal definition of variable group in MAP.'
1459 stop 'Error - illegal definition of variable group in MAP.'
1461 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1462 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1463 call readi(mapcard,'NSTEP',nstep(imap),0)
1464 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1466 & 'Illegal boundary and/or step size specification in MAP.'
1467 stop 'Illegal boundary and/or step size specification in MAP.'
1472 c----------------------------------------------------------------------------
1474 implicit real*8 (a-h,o-z)
1475 include 'DIMENSIONS'
1476 include 'COMMON.IOUNITS'
1477 include 'COMMON.GEO'
1478 include 'COMMON.CSA'
1479 include 'COMMON.BANK'
1480 include 'COMMON.CONTROL'
1482 character*620 mcmcard
1483 call card_concat(mcmcard)
1485 call readi(mcmcard,'NCONF',nconf,50)
1486 call readi(mcmcard,'NADD',nadd,0)
1487 call readi(mcmcard,'JSTART',jstart,1)
1488 call readi(mcmcard,'JEND',jend,1)
1489 call readi(mcmcard,'NSTMAX',nstmax,500000)
1490 call readi(mcmcard,'N0',n0,1)
1491 call readi(mcmcard,'N1',n1,6)
1492 call readi(mcmcard,'N2',n2,4)
1493 call readi(mcmcard,'N3',n3,0)
1494 call readi(mcmcard,'N4',n4,0)
1495 call readi(mcmcard,'N5',n5,0)
1496 call readi(mcmcard,'N6',n6,10)
1497 call readi(mcmcard,'N7',n7,0)
1498 call readi(mcmcard,'N8',n8,0)
1499 call readi(mcmcard,'N9',n9,0)
1500 call readi(mcmcard,'N14',n14,0)
1501 call readi(mcmcard,'N15',n15,0)
1502 call readi(mcmcard,'N16',n16,0)
1503 call readi(mcmcard,'N17',n17,0)
1504 call readi(mcmcard,'N18',n18,0)
1506 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1508 call readi(mcmcard,'NDIFF',ndiff,2)
1509 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1510 call readi(mcmcard,'IS1',is1,1)
1511 call readi(mcmcard,'IS2',is2,8)
1512 call readi(mcmcard,'NRAN0',nran0,4)
1513 call readi(mcmcard,'NRAN1',nran1,2)
1514 call readi(mcmcard,'IRR',irr,1)
1515 call readi(mcmcard,'NSEED',nseed,20)
1516 call readi(mcmcard,'NTOTAL',ntotal,10000)
1517 call reada(mcmcard,'CUT1',cut1,2.0d0)
1518 call reada(mcmcard,'CUT2',cut2,5.0d0)
1519 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1520 call readi(mcmcard,'ICMAX',icmax,3)
1521 call readi(mcmcard,'IRESTART',irestart,0)
1522 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1525 call reada(mcmcard,'DELE',dele,20.0d0)
1526 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1527 call readi(mcmcard,'IREF',iref,0)
1528 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1529 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1530 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1531 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1532 write (iout,*) "NCONF_IN",nconf_in
1535 c----------------------------------------------------------------------------
1536 cfmc subroutine mcmfread
1537 cfmc implicit real*8 (a-h,o-z)
1538 cfmc include 'DIMENSIONS'
1539 cfmc include 'COMMON.MCMF'
1540 cfmc include 'COMMON.IOUNITS'
1541 cfmc include 'COMMON.GEO'
1542 cfmc character*80 ucase
1543 cfmc character*620 mcmcard
1544 cfmc call card_concat(mcmcard)
1546 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1547 cfmc write(iout,*)'MAXRANT=',maxrant
1548 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1549 cfmc write(iout,*)'MAXFAM=',maxfam
1550 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1551 cfmc write(iout,*)'NNET1=',nnet1
1552 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1553 cfmc write(iout,*)'NNET2=',nnet2
1554 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1555 cfmc write(iout,*)'NNET3=',nnet3
1556 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1557 cfmc write(iout,*)'ILASTT=',ilastt
1558 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1559 cfmc write(iout,*)'MAXSTR=',maxstr
1560 cfmc maxstr_f=maxstr/maxfam
1561 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1562 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1563 cfmc write(iout,*)'NMCMF=',nmcmf
1564 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1565 cfmc write(iout,*)'IFOCUS=',ifocus
1566 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1567 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1568 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1569 cfmc write(iout,*)'INTPRT=',intprt
1570 cfmc call readi(mcmcard,'IPRT',iprt,100)
1571 cfmc write(iout,*)'IPRT=',iprt
1572 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1573 cfmc write(iout,*)'IMAXTR=',imaxtr
1574 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1575 cfmc write(iout,*)'MAXEVEN=',maxeven
1576 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1577 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1578 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1579 cfmc write(iout,*)'INIMIN=',inimin
1580 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1581 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1582 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1583 cfmc write(iout,*)'NTHREAD=',nthread
1584 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1585 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1586 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1587 cfmc write(iout,*)'MAXPERT=',maxpert
1588 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1589 cfmc write(iout,*)'IRMSD=',irmsd
1590 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1591 cfmc write(iout,*)'DENEMIN=',denemin
1592 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1593 cfmc write(iout,*)'RCUT1S=',rcut1s
1594 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1595 cfmc write(iout,*)'RCUT1E=',rcut1e
1596 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1597 cfmc write(iout,*)'RCUT2S=',rcut2s
1598 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1599 cfmc write(iout,*)'RCUT2E=',rcut2e
1600 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1601 cfmc write(iout,*)'DPERT1=',d_pert1
1602 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1603 cfmc write(iout,*)'DPERT1A=',d_pert1a
1604 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1605 cfmc write(iout,*)'DPERT2=',d_pert2
1606 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1607 cfmc write(iout,*)'DPERT2A=',d_pert2a
1608 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1609 cfmc write(iout,*)'DPERT2B=',d_pert2b
1610 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1611 cfmc write(iout,*)'DPERT2C=',d_pert2c
1612 cfmc d_pert1=deg2rad*d_pert1
1613 cfmc d_pert1a=deg2rad*d_pert1a
1614 cfmc d_pert2=deg2rad*d_pert2
1615 cfmc d_pert2a=deg2rad*d_pert2a
1616 cfmc d_pert2b=deg2rad*d_pert2b
1617 cfmc d_pert2c=deg2rad*d_pert2c
1618 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1619 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1620 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1621 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1622 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1623 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1624 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1625 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1626 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1627 cfmc write(iout,*)'RCUTINI=',rcutini
1628 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1629 cfmc write(iout,*)'GRAT=',grat
1630 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1631 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1635 c----------------------------------------------------------------------------
1637 implicit real*8 (a-h,o-z)
1638 include 'DIMENSIONS'
1639 include 'COMMON.MCM'
1640 include 'COMMON.MCE'
1641 include 'COMMON.IOUNITS'
1643 character*320 mcmcard
1644 call card_concat(mcmcard)
1645 call readi(mcmcard,'MAXACC',maxacc,100)
1646 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1647 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1648 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1649 call readi(mcmcard,'MAXREPM',maxrepm,200)
1650 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1651 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1652 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1653 call reada(mcmcard,'E_UP',e_up,5.0D0)
1654 call reada(mcmcard,'DELTE',delte,0.1D0)
1655 call readi(mcmcard,'NSWEEP',nsweep,5)
1656 call readi(mcmcard,'NSTEPH',nsteph,0)
1657 call readi(mcmcard,'NSTEPC',nstepc,0)
1658 call reada(mcmcard,'TMIN',tmin,298.0D0)
1659 call reada(mcmcard,'TMAX',tmax,298.0D0)
1660 call readi(mcmcard,'NWINDOW',nwindow,0)
1661 call readi(mcmcard,'PRINT_MC',print_mc,0)
1662 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1663 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1664 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1665 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1666 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1667 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1668 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1669 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1670 if (nwindow.gt.0) then
1671 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1673 winlen(i)=winend(i)-winstart(i)+1
1676 if (tmax.lt.tmin) tmax=tmin
1677 if (tmax.eq.tmin) then
1681 if (nstepc.gt.0 .and. nsteph.gt.0) then
1682 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1683 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1685 C Probabilities of different move types
1686 sumpro_type(0)=0.0D0
1687 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1688 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1689 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1690 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1691 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1692 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1693 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1695 print *,'i',i,' sumprotype',sumpro_type(i)
1696 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1697 print *,'i',i,' sumprotype',sumpro_type(i)
1701 c----------------------------------------------------------------------------
1702 subroutine read_minim
1703 implicit real*8 (a-h,o-z)
1704 include 'DIMENSIONS'
1705 include 'COMMON.MINIM'
1706 include 'COMMON.IOUNITS'
1708 character*320 minimcard
1709 call card_concat(minimcard)
1710 call readi(minimcard,'MAXMIN',maxmin,2000)
1711 call readi(minimcard,'MAXFUN',maxfun,5000)
1712 call readi(minimcard,'MINMIN',minmin,maxmin)
1713 call readi(minimcard,'MINFUN',minfun,maxmin)
1714 call reada(minimcard,'TOLF',tolf,1.0D-2)
1715 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1716 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1717 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1718 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1719 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1720 & 'Options in energy minimization:'
1721 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1722 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1723 & 'MinMin:',MinMin,' MinFun:',MinFun,
1724 & ' TolF:',TolF,' RTolF:',RTolF
1727 c----------------------------------------------------------------------------
1728 subroutine read_angles(kanal,*)
1729 implicit real*8 (a-h,o-z)
1730 include 'DIMENSIONS'
1731 include 'COMMON.GEO'
1732 include 'COMMON.VAR'
1733 include 'COMMON.CHAIN'
1734 include 'COMMON.IOUNITS'
1735 include 'COMMON.CONTROL'
1736 c Read angles from input
1738 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1739 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1740 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1741 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1744 c 9/7/01 avoid 180 deg valence angle
1745 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1747 theta(i)=deg2rad*theta(i)
1748 phi(i)=deg2rad*phi(i)
1749 alph(i)=deg2rad*alph(i)
1750 omeg(i)=deg2rad*omeg(i)
1755 c----------------------------------------------------------------------------
1756 subroutine reada(rekord,lancuch,wartosc,default)
1758 character*(*) rekord,lancuch
1759 double precision wartosc,default
1762 iread=index(rekord,lancuch)
1763 if (iread.eq.0) then
1767 iread=iread+ilen(lancuch)+1
1768 read (rekord(iread:),*,err=10,end=10) wartosc
1773 c----------------------------------------------------------------------------
1774 subroutine readi(rekord,lancuch,wartosc,default)
1776 character*(*) rekord,lancuch
1777 integer wartosc,default
1780 iread=index(rekord,lancuch)
1781 if (iread.eq.0) then
1785 iread=iread+ilen(lancuch)+1
1786 read (rekord(iread:),*,err=10,end=10) wartosc
1791 c----------------------------------------------------------------------------
1792 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1795 integer tablica(dim),default
1796 character*(*) rekord,lancuch
1803 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1804 if (iread.eq.0) return
1805 iread=iread+ilen(lancuch)+1
1806 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1809 c----------------------------------------------------------------------------
1810 subroutine multreada(rekord,lancuch,tablica,dim,default)
1813 double precision tablica(dim),default
1814 character*(*) rekord,lancuch
1821 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1822 if (iread.eq.0) return
1823 iread=iread+ilen(lancuch)+1
1824 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1827 c----------------------------------------------------------------------------
1828 subroutine openunits
1829 implicit real*8 (a-h,o-z)
1830 include 'DIMENSIONS'
1833 character*16 form,nodename
1836 include 'COMMON.SETUP'
1837 include 'COMMON.IOUNITS'
1839 include 'COMMON.CONTROL'
1840 integer lenpre,lenpot,ilen,lentmp
1842 character*3 out1file_text,ucase
1845 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1846 call getenv_loc("PREFIX",prefix)
1848 call getenv_loc("POT",pot)
1849 call getenv_loc("DIRTMP",tmpdir)
1850 call getenv_loc("CURDIR",curdir)
1851 call getenv_loc("OUT1FILE",out1file_text)
1852 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1853 out1file_text=ucase(out1file_text)
1854 if (out1file_text(1:1).eq."Y") then
1857 out1file=fg_rank.gt.0
1862 if (lentmp.gt.0) then
1863 write (*,'(80(1h!))')
1864 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1865 write (*,'(80(1h!))')
1866 write (*,*)"All output files will be on node /tmp directory."
1868 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1869 if (me.eq.king) then
1870 write (*,*) "The master node is ",nodename
1871 else if (fg_rank.eq.0) then
1872 write (*,*) "I am the CG slave node ",nodename
1874 write (*,*) "I am the FG slave node ",nodename
1877 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1878 lenpre = lentmp+lenpre+1
1880 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1881 C Get the names and open the input files
1882 #if defined(WINIFL) || defined(WINPGI)
1883 open(1,file=pref_orig(:ilen(pref_orig))//
1884 & '.inp',status='old',readonly,shared)
1885 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1886 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1887 C Get parameter filenames and open the parameter files.
1888 call getenv_loc('BONDPAR',bondname)
1889 open (ibond,file=bondname,status='old',readonly,shared)
1890 call getenv_loc('THETPAR',thetname)
1891 open (ithep,file=thetname,status='old',readonly,shared)
1892 call getenv_loc('ROTPAR',rotname)
1893 open (irotam,file=rotname,status='old',readonly,shared)
1894 call getenv_loc('TORPAR',torname)
1895 open (itorp,file=torname,status='old',readonly,shared)
1896 call getenv_loc('TORDPAR',tordname)
1897 open (itordp,file=tordname,status='old',readonly,shared)
1898 call getenv_loc('FOURIER',fouriername)
1899 open (ifourier,file=fouriername,status='old',readonly,shared)
1900 call getenv_loc('ELEPAR',elename)
1901 open (ielep,file=elename,status='old',readonly,shared)
1902 call getenv_loc('SIDEPAR',sidename)
1903 open (isidep,file=sidename,status='old',readonly,shared)
1904 #elif (defined CRAY) || (defined AIX)
1905 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1907 c print *,"Processor",myrank," opened file 1"
1908 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1909 c print *,"Processor",myrank," opened file 9"
1910 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1911 C Get parameter filenames and open the parameter files.
1912 call getenv_loc('BONDPAR',bondname)
1913 open (ibond,file=bondname,status='old',action='read')
1914 c print *,"Processor",myrank," opened file IBOND"
1915 call getenv_loc('THETPAR',thetname)
1916 open (ithep,file=thetname,status='old',action='read')
1917 c print *,"Processor",myrank," opened file ITHEP"
1918 call getenv_loc('ROTPAR',rotname)
1919 open (irotam,file=rotname,status='old',action='read')
1920 c print *,"Processor",myrank," opened file IROTAM"
1921 call getenv_loc('TORPAR',torname)
1922 open (itorp,file=torname,status='old',action='read')
1923 c print *,"Processor",myrank," opened file ITORP"
1924 call getenv_loc('TORDPAR',tordname)
1925 open (itordp,file=tordname,status='old',action='read')
1926 c print *,"Processor",myrank," opened file ITORDP"
1927 call getenv_loc('SCCORPAR',sccorname)
1928 open (isccor,file=sccorname,status='old',action='read')
1929 c print *,"Processor",myrank," opened file ISCCOR"
1930 call getenv_loc('FOURIER',fouriername)
1931 open (ifourier,file=fouriername,status='old',action='read')
1932 c print *,"Processor",myrank," opened file IFOURIER"
1933 call getenv_loc('ELEPAR',elename)
1934 open (ielep,file=elename,status='old',action='read')
1935 c print *,"Processor",myrank," opened file IELEP"
1936 call getenv_loc('SIDEPAR',sidename)
1937 open (isidep,file=sidename,status='old',action='read')
1938 c print *,"Processor",myrank," opened file ISIDEP"
1939 c print *,"Processor",myrank," opened parameter files"
1941 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1942 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1943 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1944 C Get parameter filenames and open the parameter files.
1945 call getenv_loc('BONDPAR',bondname)
1946 open (ibond,file=bondname,status='old')
1947 call getenv_loc('THETPAR',thetname)
1948 open (ithep,file=thetname,status='old')
1949 call getenv_loc('ROTPAR',rotname)
1950 open (irotam,file=rotname,status='old')
1951 call getenv_loc('TORPAR',torname)
1952 open (itorp,file=torname,status='old')
1953 call getenv_loc('TORDPAR',tordname)
1954 open (itordp,file=tordname,status='old')
1955 call getenv_loc('SCCORPAR',sccorname)
1956 open (isccor,file=sccorname,status='old')
1957 call getenv_loc('FOURIER',fouriername)
1958 open (ifourier,file=fouriername,status='old')
1959 call getenv_loc('ELEPAR',elename)
1960 open (ielep,file=elename,status='old')
1961 call getenv_loc('SIDEPAR',sidename)
1962 open (isidep,file=sidename,status='old')
1964 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1966 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1967 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1968 C Get parameter filenames and open the parameter files.
1969 call getenv_loc('BONDPAR',bondname)
1970 open (ibond,file=bondname,status='old',readonly)
1971 call getenv_loc('THETPAR',thetname)
1972 open (ithep,file=thetname,status='old',readonly)
1973 call getenv_loc('ROTPAR',rotname)
1974 open (irotam,file=rotname,status='old',readonly)
1975 call getenv_loc('TORPAR',torname)
1976 open (itorp,file=torname,status='old',readonly)
1977 call getenv_loc('TORDPAR',tordname)
1978 open (itordp,file=tordname,status='old',readonly)
1979 call getenv_loc('SCCORPAR',sccorname)
1980 open (isccor,file=sccorname,status='old',readonly)
1982 call getenv_loc('THETPARPDB',thetname_pdb)
1983 print *,"thetname_pdb ",thetname_pdb
1984 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1985 print *,ithep_pdb," opened"
1987 call getenv_loc('FOURIER',fouriername)
1988 open (ifourier,file=fouriername,status='old',readonly)
1989 call getenv_loc('ELEPAR',elename)
1990 open (ielep,file=elename,status='old',readonly)
1991 call getenv_loc('SIDEPAR',sidename)
1992 open (isidep,file=sidename,status='old',readonly)
1994 call getenv_loc('ROTPARPDB',rotname_pdb)
1995 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2000 C 8/9/01 In the newest version SCp interaction constants are read from a file
2001 C Use -DOLDSCP to use hard-coded constants instead.
2003 call getenv_loc('SCPPAR',scpname)
2004 #if defined(WINIFL) || defined(WINPGI)
2005 open (iscpp,file=scpname,status='old',readonly,shared)
2006 #elif (defined CRAY) || (defined AIX)
2007 open (iscpp,file=scpname,status='old',action='read')
2009 open (iscpp,file=scpname,status='old')
2011 open (iscpp,file=scpname,status='old',readonly)
2014 call getenv_loc('PATTERN',patname)
2015 #if defined(WINIFL) || defined(WINPGI)
2016 open (icbase,file=patname,status='old',readonly,shared)
2017 #elif (defined CRAY) || (defined AIX)
2018 open (icbase,file=patname,status='old',action='read')
2020 open (icbase,file=patname,status='old')
2022 open (icbase,file=patname,status='old',readonly)
2025 C Open output file only for CG processes
2026 c print *,"Processor",myrank," fg_rank",fg_rank
2027 if (fg_rank.eq.0) then
2029 if (nodes.eq.1) then
2032 npos = dlog10(dfloat(nodes-1))+1
2034 if (npos.lt.3) npos=3
2035 write (liczba,'(i1)') npos
2036 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2038 write (liczba,form) me
2039 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2040 & liczba(:ilen(liczba))
2041 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2043 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2045 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2046 & liczba(:ilen(liczba))//'.mol2'
2047 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2048 & liczba(:ilen(liczba))//'.stat'
2050 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2051 & //liczba(:ilen(liczba))//'.stat')
2052 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2055 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2056 & liczba(:ilen(liczba))//'.const'
2061 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2062 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2063 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2064 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2065 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2067 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2069 rest2name=prefix(:ilen(prefix))//'.rst'
2071 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2074 #if defined(AIX) || defined(PGI)
2075 if (me.eq.king .or. .not. out1file)
2076 & open(iout,file=outname,status='unknown')
2078 if (fg_rank.gt.0) then
2079 write (liczba,'(i3.3)') myrank/nfgtasks
2080 write (ll,'(bz,i3.3)') fg_rank
2081 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2086 open(igeom,file=intname,status='unknown',position='append')
2087 open(ipdb,file=pdbname,status='unknown')
2088 open(imol2,file=mol2name,status='unknown')
2089 open(istat,file=statname,status='unknown',position='append')
2091 c1out open(iout,file=outname,status='unknown')
2094 if (me.eq.king .or. .not.out1file)
2095 & open(iout,file=outname,status='unknown')
2097 if (fg_rank.gt.0) then
2098 write (liczba,'(i3.3)') myrank/nfgtasks
2099 write (ll,'(bz,i3.3)') fg_rank
2100 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2105 open(igeom,file=intname,status='unknown',access='append')
2106 open(ipdb,file=pdbname,status='unknown')
2107 open(imol2,file=mol2name,status='unknown')
2108 open(istat,file=statname,status='unknown',access='append')
2110 c1out open(iout,file=outname,status='unknown')
2113 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2114 csa_seed=prefix(:lenpre)//'.CSA.seed'
2115 csa_history=prefix(:lenpre)//'.CSA.history'
2116 csa_bank=prefix(:lenpre)//'.CSA.bank'
2117 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2118 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2119 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2120 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2121 csa_int=prefix(:lenpre)//'.int'
2122 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2123 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2124 csa_in=prefix(:lenpre)//'.CSA.in'
2125 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2128 write (iout,'(80(1h-))')
2129 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2130 write (iout,'(80(1h-))')
2131 write (iout,*) "Input file : ",
2132 & pref_orig(:ilen(pref_orig))//'.inp'
2133 write (iout,*) "Output file : ",
2134 & outname(:ilen(outname))
2136 write (iout,*) "Sidechain potential file : ",
2137 & sidename(:ilen(sidename))
2139 write (iout,*) "SCp potential file : ",
2140 & scpname(:ilen(scpname))
2142 write (iout,*) "Electrostatic potential file : ",
2143 & elename(:ilen(elename))
2144 write (iout,*) "Cumulant coefficient file : ",
2145 & fouriername(:ilen(fouriername))
2146 write (iout,*) "Torsional parameter file : ",
2147 & torname(:ilen(torname))
2148 write (iout,*) "Double torsional parameter file : ",
2149 & tordname(:ilen(tordname))
2150 write (iout,*) "SCCOR parameter file : ",
2151 & sccorname(:ilen(sccorname))
2152 write (iout,*) "Bond & inertia constant file : ",
2153 & bondname(:ilen(bondname))
2154 write (iout,*) "Bending parameter file : ",
2155 & thetname(:ilen(thetname))
2156 write (iout,*) "Rotamer parameter file : ",
2157 & rotname(:ilen(rotname))
2158 write (iout,*) "Thetpdb parameter file : ",
2159 & thetname_pdb(:ilen(thetname_pdb))
2160 write (iout,*) "Threading database : ",
2161 & patname(:ilen(patname))
2163 &write (iout,*)" DIRTMP : ",
2165 write (iout,'(80(1h-))')
2169 c----------------------------------------------------------------------------
2170 subroutine card_concat(card)
2171 implicit real*8 (a-h,o-z)
2172 include 'DIMENSIONS'
2173 include 'COMMON.IOUNITS'
2175 character*80 karta,ucase
2177 read (inp,'(a)') karta
2180 do while (karta(80:80).eq.'&')
2181 card=card(:ilen(card)+1)//karta(:79)
2182 read (inp,'(a)') karta
2185 card=card(:ilen(card)+1)//karta
2188 c----------------------------------------------------------------------------------
2190 implicit real*8 (a-h,o-z)
2191 include 'DIMENSIONS'
2192 include 'COMMON.CHAIN'
2193 include 'COMMON.IOUNITS'
2195 open(irest2,file=rest2name,status='unknown')
2196 read(irest2,*) totT,EK,potE,totE,t_bath
2198 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2201 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2204 read (irest2,*) iset
2209 c---------------------------------------------------------------------------------
2210 subroutine read_fragments
2211 implicit real*8 (a-h,o-z)
2212 include 'DIMENSIONS'
2216 include 'COMMON.SETUP'
2217 include 'COMMON.CHAIN'
2218 include 'COMMON.IOUNITS'
2220 include 'COMMON.CONTROL'
2221 read(inp,*) nset,nfrag,npair,nfrag_back
2222 if(me.eq.king.or..not.out1file)
2223 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2224 & " nfrag_back",nfrag_back
2226 read(inp,*) mset(iset)
2228 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2230 if(me.eq.king.or..not.out1file)
2231 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2232 & ifrag(2,i,iset), qinfrag(i,iset)
2235 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2237 if(me.eq.king.or..not.out1file)
2238 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2239 & ipair(2,i,iset), qinpair(i,iset)
2242 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2243 & wfrag_back(3,i,iset),
2244 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2245 if(me.eq.king.or..not.out1file)
2246 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2247 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2252 c-------------------------------------------------------------------------------
2253 subroutine read_dist_constr
2254 implicit real*8 (a-h,o-z)
2255 include 'DIMENSIONS'
2259 include 'COMMON.SETUP'
2260 include 'COMMON.CONTROL'
2261 include 'COMMON.CHAIN'
2262 include 'COMMON.IOUNITS'
2263 include 'COMMON.SBRIDGE'
2264 integer ifrag_(2,100),ipair_(2,100)
2265 double precision wfrag_(100),wpair_(100)
2266 character*500 controlcard
2268 write (iout,*) "Calling read_dist_constr"
2269 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2271 call card_concat(controlcard)
2272 call readi(controlcard,"NFRAG",nfrag_,0)
2273 call readi(controlcard,"NPAIR",npair_,0)
2274 call readi(controlcard,"NDIST",ndist_,0)
2275 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2276 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2277 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2278 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2279 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2280 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2281 c write (iout,*) "IFRAG"
2283 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2285 c write (iout,*) "IPAIR"
2287 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2291 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2292 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2293 & ifrag_(2,i)=nstart_sup+nsup-1
2294 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2296 if (wfrag_(i).gt.0.0d0) then
2297 do j=ifrag_(1,i),ifrag_(2,i)-1
2298 do k=j+1,ifrag_(2,i)
2299 c write (iout,*) "j",j," k",k
2301 if (constr_dist.eq.1) then
2306 forcon(nhpb)=wfrag_(i)
2307 else if (constr_dist.eq.2) then
2308 if (ddjk.le.dist_cut) then
2313 forcon(nhpb)=wfrag_(i)
2320 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2323 if (.not.out1file .or. me.eq.king)
2324 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2325 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2327 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2328 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2335 if (wpair_(i).gt.0.0d0) then
2343 do j=ifrag_(1,ii),ifrag_(2,ii)
2344 do k=ifrag_(1,jj),ifrag_(2,jj)
2348 forcon(nhpb)=wpair_(i)
2349 dhpb(nhpb)=dist(j,k)
2351 if (.not.out1file .or. me.eq.king)
2352 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2353 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2355 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2356 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2364 if (constr_dist.eq.11) then
2365 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2366 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2367 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2370 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2371 & ibecarb(i),forcon(nhpb+1)
2373 if (forcon(nhpb+1).gt.0.0d0) then
2375 if (ibecarb(i).gt.0) then
2376 ihpb(i)=ihpb(i)+nres
2377 jhpb(i)=jhpb(i)+nres
2379 if (dhpb(nhpb).eq.0.0d0)
2380 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2382 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2383 C if (forcon(nhpb+1).gt.0.0d0) then
2385 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2387 if (.not.out1file .or. me.eq.king)
2388 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2389 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2391 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2392 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2399 c-------------------------------------------------------------------------------
2401 subroutine flush(iu)
2406 subroutine flush(iu)
2411 c------------------------------------------------------------------------------
2412 subroutine copy_to_tmp(source)
2413 include "DIMENSIONS"
2414 include "COMMON.IOUNITS"
2415 character*(*) source
2416 character* 256 tmpfile
2420 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2421 inquire(file=tmpfile,exist=ex)
2423 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2424 & " to temporary directory..."
2425 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2426 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2430 c------------------------------------------------------------------------------
2431 subroutine move_from_tmp(source)
2432 include "DIMENSIONS"
2433 include "COMMON.IOUNITS"
2434 character*(*) source
2437 write (*,*) "Moving ",source(:ilen(source)),
2438 & " from temporary directory to working directory"
2439 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2440 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2443 c------------------------------------------------------------------------------
2444 subroutine random_init(seed)
2446 C Initialize random number generator
2448 implicit real*8 (a-h,o-z)
2449 include 'DIMENSIONS'
2452 logical OKRandom, prng_restart
2454 integer iseed_array(4)
2456 include 'COMMON.IOUNITS'
2457 include 'COMMON.TIME1'
2458 include 'COMMON.THREAD'
2459 include 'COMMON.SBRIDGE'
2460 include 'COMMON.CONTROL'
2461 include 'COMMON.MCM'
2462 include 'COMMON.MAP'
2463 include 'COMMON.HEADER'
2464 include 'COMMON.CSA'
2465 include 'COMMON.CHAIN'
2466 include 'COMMON.MUCA'
2468 include 'COMMON.FFIELD'
2469 include 'COMMON.SETUP'
2470 iseed=-dint(dabs(seed))
2471 if (iseed.eq.0) then
2472 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2473 & 'Random seed undefined. The program will stop.'
2474 write (*,'(/80(1h*)/20x,a/80(1h*))')
2475 & 'Random seed undefined. The program will stop.'
2477 call mpi_finalize(mpi_comm_world,ierr)
2479 stop 'Bad random seed.'
2482 if (fg_rank.eq.0) then
2486 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2487 OKRandom = prng_restart(me,iseed)
2490 tmp=65536.0d0**(4-i)
2491 iseed_array(i) = dint(seed/tmp)
2492 seed=seed-iseed_array(i)*tmp
2495 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2496 & (iseed_array(i),i=1,4)
2497 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2498 & (iseed_array(i),i=1,4)
2499 OKRandom = prng_restart(me,iseed_array)
2502 c r1 = prng_next(me)
2503 r1=ran_number(0.0D0,1.0D0)
2505 & write (iout,*) 'ran_num',r1
2506 if (r1.lt.0.0d0) OKRandom=.false.
2508 if (.not.OKRandom) then
2509 write (iout,*) 'PRNG IS NOT WORKING!!!'
2510 print *,'PRNG IS NOT WORKING!!!'
2513 call mpi_abort(mpi_comm_world,error_msg,ierr)
2516 write (iout,*) 'too many processors for parallel prng'
2517 write (*,*) 'too many processors for parallel prng'
2525 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)