2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 call readi(controlcard,'SYM',symetr,1)
102 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
103 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
104 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
105 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
106 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
107 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
108 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
109 call reada(controlcard,'DRMS',drms,0.1D0)
110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
111 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
112 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
113 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
114 write (iout,'(a,f10.1)')'DRMS = ',drms
115 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
116 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
118 call readi(controlcard,'NZ_START',nz_start,0)
119 call readi(controlcard,'NZ_END',nz_end,0)
120 call readi(controlcard,'IZ_SC',iz_sc,0)
122 safety = 60.0d0*safety
125 call reada(controlcard,"T_BATH",t_bath,300.0d0)
126 minim=(index(controlcard,'MINIMIZE').gt.0)
127 dccart=(index(controlcard,'CART').gt.0)
128 overlapsc=(index(controlcard,'OVERLAP').gt.0)
129 overlapsc=.not.overlapsc
130 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
131 searchsc=.not.searchsc
132 sideadd=(index(controlcard,'SIDEADD').gt.0)
133 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
134 outpdb=(index(controlcard,'PDBOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 AFMlog=(index(controlcard,'AFM'))
141 selfguide=(index(controlcard,'SELFGUIDE'))
142 print *,'AFMlog',AFMlog,selfguide,"KUPA"
143 call readi(controlcard,'IPRINT',iprint,0)
144 call readi(controlcard,'MAXGEN',maxgen,10000)
145 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
146 call readi(controlcard,"KDIAG",kdiag,0)
147 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
148 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
149 & write (iout,*) "RESCALE_MODE",rescale_mode
150 split_ene=index(controlcard,'SPLIT_ENE').gt.0
151 if (index(controlcard,'REGULAR').gt.0.0D0) then
152 call reada(controlcard,'WEIDIS',weidis,0.1D0)
156 if (index(controlcard,'CHECKGRAD').gt.0) then
158 if (index(controlcard,'CART').gt.0) then
160 elseif (index(controlcard,'CARINT').gt.0) then
165 elseif (index(controlcard,'THREAD').gt.0) then
167 call readi(controlcard,'THREAD',nthread,0)
168 if (nthread.gt.0) then
169 call reada(controlcard,'WEIDIS',weidis,0.1D0)
172 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
173 stop 'Error termination in Read_Control.'
175 else if (index(controlcard,'MCMA').gt.0) then
177 else if (index(controlcard,'MCEE').gt.0) then
179 else if (index(controlcard,'MULTCONF').gt.0) then
181 else if (index(controlcard,'MAP').gt.0) then
183 call readi(controlcard,'MAP',nmap,0)
184 else if (index(controlcard,'CSA').gt.0) then
186 crc else if (index(controlcard,'ZSCORE').gt.0) then
188 crc ZSCORE is rm from UNRES, modecalc=9 is available
191 cfcm else if (index(controlcard,'MCMF').gt.0) then
193 else if (index(controlcard,'SOFTREG').gt.0) then
195 else if (index(controlcard,'CHECK_BOND').gt.0) then
197 else if (index(controlcard,'TEST').gt.0) then
199 else if (index(controlcard,'MD').gt.0) then
201 else if (index(controlcard,'RE ').gt.0) then
205 lmuca=index(controlcard,'MUCA').gt.0
206 call readi(controlcard,'MUCADYN',mucadyn,0)
207 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
208 if (lmuca .and. (me.eq.king .or. .not.out1file ))
210 write (iout,*) 'MUCADYN=',mucadyn
211 write (iout,*) 'MUCASMOOTH=',muca_smooth
214 iscode=index(controlcard,'ONE_LETTER')
215 indphi=index(controlcard,'PHI')
216 indback=index(controlcard,'BACK')
217 iranconf=index(controlcard,'RAND_CONF')
218 i2ndstr=index(controlcard,'USE_SEC_PRED')
219 gradout=index(controlcard,'GRADOUT').gt.0
220 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
221 C DISTCHAINMAX become obsolete for periodic boundry condition
222 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
223 C Reading the dimensions of box in x,y,z coordinates
224 call reada(controlcard,'BOXX',boxxsize,100.0d0)
225 call reada(controlcard,'BOXY',boxysize,100.0d0)
226 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
227 c Cutoff range for interactions
228 call reada(controlcard,"R_CUT",r_cut,15.0d0)
229 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
230 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
231 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
232 if (lipthick.gt.0.0d0) then
233 bordliptop=(boxzsize+lipthick)/2.0
234 bordlipbot=bordliptop-lipthick
236 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
237 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
238 buflipbot=bordlipbot+lipbufthick
239 bufliptop=bordliptop-lipbufthick
240 if ((lipbufthick*2.0d0).gt.lipthick)
241 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
243 write(iout,*) "bordliptop=",bordliptop
244 write(iout,*) "bordlipbot=",bordlipbot
245 write(iout,*) "bufliptop=",bufliptop
246 write(iout,*) "buflipbot=",buflipbot
249 if (me.eq.king .or. .not.out1file )
250 & write (iout,*) "DISTCHAINMAX",distchainmax
252 if(me.eq.king.or..not.out1file)
253 & write (iout,'(2a)') diagmeth(kdiag),
254 & ' routine used to diagonalize matrices.'
257 c--------------------------------------------------------------------------
258 subroutine read_REMDpar
262 implicit real*8 (a-h,o-z)
264 include 'COMMON.IOUNITS'
265 include 'COMMON.TIME1'
268 include 'COMMON.LANGEVIN'
270 include 'COMMON.LANGEVIN.lang0'
272 include 'COMMON.INTERACT'
273 include 'COMMON.NAMES'
275 include 'COMMON.REMD'
276 include 'COMMON.CONTROL'
277 include 'COMMON.SETUP'
279 character*320 controlcard
280 character*3200 controlcard1
281 integer iremd_m_total
283 if(me.eq.king.or..not.out1file)
284 & write (iout,*) "REMD setup"
286 call card_concat(controlcard)
287 call readi(controlcard,"NREP",nrep,3)
288 call readi(controlcard,"NSTEX",nstex,1000)
289 call reada(controlcard,"RETMIN",retmin,10.0d0)
290 call reada(controlcard,"RETMAX",retmax,1000.0d0)
291 mremdsync=(index(controlcard,'SYNC').gt.0)
292 call readi(controlcard,"NSYN",i_sync_step,100)
293 restart1file=(index(controlcard,'REST1FILE').gt.0)
294 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
295 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
296 if(max_cache_traj_use.gt.max_cache_traj)
297 & max_cache_traj_use=max_cache_traj
298 if(me.eq.king.or..not.out1file) then
299 cd if (traj1file) then
300 crc caching is in testing - NTWX is not ignored
301 cd write (iout,*) "NTWX value is ignored"
302 cd write (iout,*) " trajectory is stored to one file by master"
303 cd write (iout,*) " before exchange at NSTEX intervals"
305 write (iout,*) "NREP= ",nrep
306 write (iout,*) "NSTEX= ",nstex
307 write (iout,*) "SYNC= ",mremdsync
308 write (iout,*) "NSYN= ",i_sync_step
309 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
312 if (index(controlcard,'TLIST').gt.0) then
314 call card_concat(controlcard1)
315 read(controlcard1,*) (remd_t(i),i=1,nrep)
316 if(me.eq.king.or..not.out1file)
317 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
320 if (index(controlcard,'MLIST').gt.0) then
322 call card_concat(controlcard1)
323 read(controlcard1,*) (remd_m(i),i=1,nrep)
324 if(me.eq.king.or..not.out1file) then
325 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
328 iremd_m_total=iremd_m_total+remd_m(i)
330 write (iout,*) 'Total number of replicas ',iremd_m_total
333 if(me.eq.king.or..not.out1file)
334 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
337 c--------------------------------------------------------------------------
338 subroutine read_MDpar
342 implicit real*8 (a-h,o-z)
344 include 'COMMON.IOUNITS'
345 include 'COMMON.TIME1'
348 include 'COMMON.LANGEVIN'
350 include 'COMMON.LANGEVIN.lang0'
352 include 'COMMON.INTERACT'
353 include 'COMMON.NAMES'
355 include 'COMMON.SETUP'
356 include 'COMMON.CONTROL'
357 include 'COMMON.SPLITELE'
359 character*320 controlcard
361 call card_concat(controlcard)
362 call readi(controlcard,"NSTEP",n_timestep,1000000)
363 call readi(controlcard,"NTWE",ntwe,100)
364 call readi(controlcard,"NTWX",ntwx,1000)
365 call reada(controlcard,"DT",d_time,1.0d-1)
366 call reada(controlcard,"DVMAX",dvmax,2.0d1)
367 call reada(controlcard,"DAMAX",damax,1.0d1)
368 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
369 call readi(controlcard,"LANG",lang,0)
370 RESPA = index(controlcard,"RESPA") .gt. 0
371 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
372 ntime_split0=ntime_split
373 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
374 ntime_split0=ntime_split
375 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
376 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
377 rest = index(controlcard,"REST").gt.0
378 tbf = index(controlcard,"TBF").gt.0
379 usampl = index(controlcard,"USAMPL").gt.0
381 mdpdb = index(controlcard,"MDPDB").gt.0
382 call reada(controlcard,"T_BATH",t_bath,300.0d0)
383 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
384 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
385 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
386 if (count_reset_moment.eq.0) count_reset_moment=1000000000
387 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
388 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
389 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
390 if (count_reset_vel.eq.0) count_reset_vel=1000000000
391 large = index(controlcard,"LARGE").gt.0
392 print_compon = index(controlcard,"PRINT_COMPON").gt.0
393 rattle = index(controlcard,"RATTLE").gt.0
394 preminim = index(controlcard,"PREMINIM").gt.0
396 dccart=(index(controlcard,'CART').gt.0)
399 c if performing umbrella sampling, fragments constrained are read from the fragment file
405 if(me.eq.king.or..not.out1file) then
407 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
409 write (iout,'(a)') "The units are:"
410 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
411 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
412 & " acceleration: angstrom/(48.9 fs)**2"
413 write (iout,'(a)') "energy: kcal/mol, temperature: K"
415 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
416 write (iout,'(a60,f10.5,a)')
417 & "Initial time step of numerical integration:",d_time,
419 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
421 write (iout,'(2a,i4,a)')
422 & "A-MTS algorithm used; initial time step for fast-varying",
423 & " short-range forces split into",ntime_split," steps."
424 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
425 & r_cut," lambda",rlamb
427 write (iout,'(2a,f10.5)')
428 & "Maximum acceleration threshold to reduce the time step",
429 & "/increase split number:",damax
430 write (iout,'(2a,f10.5)')
431 & "Maximum predicted energy drift to reduce the timestep",
432 & "/increase split number:",edriftmax
433 write (iout,'(a60,f10.5)')
434 & "Maximum velocity threshold to reduce velocities:",dvmax
435 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
436 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
437 if (rattle) write (iout,'(a60)')
438 & "Rattle algorithm used to constrain the virtual bonds"
439 if (preminim .or. iranconf.gt.0) then
441 & "Initial structure will be energy-minimized"
446 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
447 call reada(controlcard,"RWAT",rwat,1.4d0)
448 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
449 surfarea=index(controlcard,"SURFAREA").gt.0
450 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
451 if(me.eq.king.or..not.out1file)then
452 write (iout,'(/a,$)') "Langevin dynamics calculation"
455 & " with direct integration of Langevin equations"
456 else if (lang.eq.2) then
457 write (iout,'(a/)') " with TINKER stochasic MD integrator"
458 else if (lang.eq.3) then
459 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
460 else if (lang.eq.4) then
461 write (iout,'(a/)') " in overdamped mode"
463 write (iout,'(//a,i5)')
464 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
467 write (iout,'(a60,f10.5)') "Temperature:",t_bath
468 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
469 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
470 write (iout,'(a60,f10.5)')
471 & "Scaling factor of the friction forces:",scal_fric
472 if (surfarea) write (iout,'(2a,i10,a)')
473 & "Friction coefficients will be scaled by solvent-accessible",
474 & " surface area every",reset_fricmat," steps."
476 c Calculate friction coefficients and bounds of stochastic forces
477 eta=6*pi*cPoise*etawat
478 if(me.eq.king.or..not.out1file)
479 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
481 gamp=scal_fric*(pstok+rwat)*eta
482 stdfp=dsqrt(2*Rb*t_bath/d_time)
484 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
485 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
487 if(me.eq.king.or..not.out1file)then
488 write (iout,'(/2a/)')
489 & "Radii of site types and friction coefficients and std's of",
490 & " stochastic forces of fully exposed sites"
491 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
493 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
494 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
498 if(me.eq.king.or..not.out1file)then
499 write (iout,'(a)') "Berendsen bath calculation"
500 write (iout,'(a60,f10.5)') "Temperature:",t_bath
501 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
503 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
504 & count_reset_moment," steps"
506 & write (iout,'(a,i10,a)')
507 & "Velocities will be reset at random every",count_reset_vel,
511 if(me.eq.king.or..not.out1file)
512 & write (iout,'(a31)') "Microcanonical mode calculation"
514 if(me.eq.king.or..not.out1file)then
515 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
517 write(iout,*) "MD running with constraints."
518 write(iout,*) "Equilibration time ", eq_time, " mtus."
519 write(iout,*) "Constraining ", nfrag," fragments."
520 write(iout,*) "Length of each fragment, weight and q0:"
522 write (iout,*) "Set of restraints #",iset
524 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
525 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
527 write(iout,*) "constraints between ", npair, "fragments."
528 write(iout,*) "constraint pairs, weights and q0:"
530 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
531 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
533 write(iout,*) "angle constraints within ", nfrag_back,
534 & "backbone fragments."
535 write(iout,*) "fragment, weights:"
537 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
538 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
539 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
542 iset=mod(kolor,nset)+1
545 if(me.eq.king.or..not.out1file)
546 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
549 c------------------------------------------------------------------------------
552 C Read molecular data.
554 implicit real*8 (a-h,o-z)
560 include 'COMMON.IOUNITS'
563 include 'COMMON.INTERACT'
564 include 'COMMON.LOCAL'
565 include 'COMMON.NAMES'
566 include 'COMMON.CHAIN'
567 include 'COMMON.FFIELD'
568 include 'COMMON.SBRIDGE'
569 include 'COMMON.HEADER'
570 include 'COMMON.CONTROL'
571 include 'COMMON.DBASE'
572 include 'COMMON.THREAD'
573 include 'COMMON.CONTACTS'
574 include 'COMMON.TORCNSTR'
575 include 'COMMON.TIME1'
576 include 'COMMON.BOUNDS'
578 include 'COMMON.SETUP'
579 character*4 sequence(maxres)
581 double precision x(maxvar)
582 character*256 pdbfile
583 character*320 weightcard
584 character*80 weightcard_t,ucase
585 dimension itype_pdb(maxres)
586 common /pizda/ itype_pdb
587 logical seq_comp,fail
588 double precision energia(0:n_ene)
594 C Read weights of the subsequent energy terms.
595 call card_concat(weightcard)
596 call reada(weightcard,'WLONG',wlong,1.0D0)
597 call reada(weightcard,'WSC',wsc,wlong)
598 call reada(weightcard,'WSCP',wscp,wlong)
599 call reada(weightcard,'WELEC',welec,1.0D0)
600 call reada(weightcard,'WVDWPP',wvdwpp,welec)
601 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
602 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
603 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
604 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
605 call reada(weightcard,'WTURN3',wturn3,1.0D0)
606 call reada(weightcard,'WTURN4',wturn4,1.0D0)
607 call reada(weightcard,'WTURN6',wturn6,1.0D0)
608 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
609 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
610 call reada(weightcard,'WBOND',wbond,1.0D0)
611 call reada(weightcard,'WTOR',wtor,1.0D0)
612 call reada(weightcard,'WTORD',wtor_d,1.0D0)
613 call reada(weightcard,'WANG',wang,1.0D0)
614 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
615 call reada(weightcard,'SCAL14',scal14,0.4D0)
616 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
617 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
618 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
619 call reada(weightcard,'TEMP0',temp0,300.0d0)
620 call reada(weightcard,'WLT',wliptran,0.0D0)
621 if (index(weightcard,'SOFT').gt.0) ipot=6
622 C 12/1/95 Added weight for the multi-body term WCORR
623 call reada(weightcard,'WCORRH',wcorr,1.0D0)
624 if (wcorr4.gt.0.0d0) wcorr=wcorr4
644 if(me.eq.king.or..not.out1file)
645 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
646 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
648 10 format (/'Energy-term weights (unscaled):'//
649 & 'WSCC= ',f10.6,' (SC-SC)'/
650 & 'WSCP= ',f10.6,' (SC-p)'/
651 & 'WELEC= ',f10.6,' (p-p electr)'/
652 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
653 & 'WBOND= ',f10.6,' (stretching)'/
654 & 'WANG= ',f10.6,' (bending)'/
655 & 'WSCLOC= ',f10.6,' (SC local)'/
656 & 'WTOR= ',f10.6,' (torsional)'/
657 & 'WTORD= ',f10.6,' (double torsional)'/
658 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
659 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
660 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
661 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
662 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
663 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
664 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
665 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
666 & 'WTURN6= ',f10.6,' (turns, 6th order)')
667 if(me.eq.king.or..not.out1file)then
668 if (wcorr4.gt.0.0d0) then
669 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
670 & 'between contact pairs of peptide groups'
671 write (iout,'(2(a,f5.3/))')
672 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
673 & 'Range of quenching the correlation terms:',2*delt_corr
674 else if (wcorr.gt.0.0d0) then
675 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
676 & 'between contact pairs of peptide groups'
678 write (iout,'(a,f8.3)')
679 & 'Scaling factor of 1,4 SC-p interactions:',scal14
680 write (iout,'(a,f8.3)')
681 & 'General scaling factor of SC-p interactions:',scalscp
683 r0_corr=cutoff_corr-delt_corr
685 aad(i,1)=scalscp*aad(i,1)
686 aad(i,2)=scalscp*aad(i,2)
687 bad(i,1)=scalscp*bad(i,1)
688 bad(i,2)=scalscp*bad(i,2)
690 call rescale_weights(t_bath)
691 if(me.eq.king.or..not.out1file)
692 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
693 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
695 22 format (/'Energy-term weights (scaled):'//
696 & 'WSCC= ',f10.6,' (SC-SC)'/
697 & 'WSCP= ',f10.6,' (SC-p)'/
698 & 'WELEC= ',f10.6,' (p-p electr)'/
699 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
700 & 'WBOND= ',f10.6,' (stretching)'/
701 & 'WANG= ',f10.6,' (bending)'/
702 & 'WSCLOC= ',f10.6,' (SC local)'/
703 & 'WTOR= ',f10.6,' (torsional)'/
704 & 'WTORD= ',f10.6,' (double torsional)'/
705 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
706 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
707 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
708 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
709 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
710 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
711 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
712 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
713 & 'WTURN6= ',f10.6,' (turns, 6th order)')
714 if(me.eq.king.or..not.out1file)
715 & write (iout,*) "Reference temperature for weights calculation:",
717 call reada(weightcard,"D0CM",d0cm,3.78d0)
718 call reada(weightcard,"AKCM",akcm,15.1d0)
719 call reada(weightcard,"AKTH",akth,11.0d0)
720 call reada(weightcard,"AKCT",akct,12.0d0)
721 call reada(weightcard,"V1SS",v1ss,-1.08d0)
722 call reada(weightcard,"V2SS",v2ss,7.61d0)
723 call reada(weightcard,"V3SS",v3ss,13.7d0)
724 call reada(weightcard,"EBR",ebr,-5.50D0)
725 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
727 dyn_ss_mask(i)=.false.
731 dyn_ssbond_ij(i,j)=1.0d300
734 call reada(weightcard,"HT",Ht,0.0D0)
736 ss_depth=ebr/wsc-0.25*eps(1,1)
737 Ht=Ht/wsc-0.25*eps(1,1)
738 akcm=akcm*wstrain/wsc
739 akth=akth*wstrain/wsc
740 akct=akct*wstrain/wsc
741 v1ss=v1ss*wstrain/wsc
742 v2ss=v2ss*wstrain/wsc
743 v3ss=v3ss*wstrain/wsc
745 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
748 if(me.eq.king.or..not.out1file) then
749 write (iout,*) "Parameters of the SS-bond potential:"
750 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
752 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
753 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
754 write (iout,*)" HT",Ht
755 print *,'indpdb=',indpdb,' pdbref=',pdbref
757 if (indpdb.gt.0 .or. pdbref) then
758 read(inp,'(a)') pdbfile
759 if(me.eq.king.or..not.out1file)
760 & write (iout,'(2a)') 'PDB data will be read from file ',
761 & pdbfile(:ilen(pdbfile))
762 open(ipdbin,file=pdbfile,status='old',err=33)
764 33 write (iout,'(a)') 'Error opening PDB file.'
767 c print *,'Begin reading pdb data'
769 c print *,'Finished reading pdb data'
770 if(me.eq.king.or..not.out1file)
771 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
772 & ' nstart_sup=',nstart_sup
774 itype_pdb(i)=itype(i)
778 nct=nstart_sup+nsup-1
779 call contact(.false.,ncont_ref,icont_ref,co)
782 if(me.eq.king.or..not.out1file)
783 & write(iout,*)'Adding sidechains'
787 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
790 do while (fail.and.nsi.le.maxsi)
791 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
794 if(fail) write(iout,*)'Adding sidechain failed for res ',
795 & i,' after ',nsi,' trials'
800 if (indpdb.eq.0) then
801 C Read sequence if not taken from the pdb file.
803 c print *,'nres=',nres
804 if (iscode.gt.0) then
805 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
807 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
809 C Convert sequence to numeric code
811 itype(i)=rescode(i,sequence(i),iscode)
813 C Assign initial virtual bond lengths
819 vbld(i+nres)=dsc(iabs(itype(i)))
820 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
821 c write (iout,*) "i",i," itype",itype(i),
822 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
826 c print '(20i4)',(itype(i),i=1,nres)
829 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
831 if (itype(i).eq.ntyp1) then
835 else if (iabs(itype(i+1)).ne.20) then
837 else if (iabs(itype(i)).ne.20) then
844 if(me.eq.king.or..not.out1file)then
845 write (iout,*) "ITEL"
847 write (iout,*) i,itype(i),itel(i)
849 print *,'Call Read_Bridge.'
852 C 8/13/98 Set limits to generating the dihedral angles
857 read (inp,*) ndih_constr
858 if (ndih_constr.gt.0) then
860 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
861 if(me.eq.king.or..not.out1file)then
863 & 'There are',ndih_constr,' constraints on phi angles.'
865 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
869 phi0(i)=deg2rad*phi0(i)
870 drange(i)=deg2rad*drange(i)
872 if(me.eq.king.or..not.out1file)
873 & write (iout,*) 'FTORS',ftors
876 phibound(1,ii) = phi0(i)-drange(i)
877 phibound(2,ii) = phi0(i)+drange(i)
884 write (iout,'(a)') 'Boundaries in phi angle sampling:'
886 write (iout,'(a3,i5,2f10.1)')
887 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
893 cd print *,'NNT=',NNT,' NCT=',NCT
894 if (itype(1).eq.ntyp1) nnt=2
895 if (itype(nres).eq.ntyp1) nct=nct-1
897 if(me.eq.king.or..not.out1file)
898 & write (iout,'(a,i3)') 'nsup=',nsup
900 if (nsup.le.(nct-nnt+1)) then
901 do i=0,nct-nnt+1-nsup
902 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
908 & 'Error - sequences to be superposed do not match.'
911 do i=0,nsup-(nct-nnt+1)
912 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
914 nstart_sup=nstart_sup+i
920 & 'Error - sequences to be superposed do not match.'
923 if (nsup.eq.0) nsup=nct-nnt
924 if (nstart_sup.eq.0) nstart_sup=nnt
925 if (nstart_seq.eq.0) nstart_seq=nnt
926 if(me.eq.king.or..not.out1file)
927 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
928 & ' nstart_seq=',nstart_seq
930 c--- Zscore rms -------
931 if (nz_start.eq.0) nz_start=nnt
932 if (nz_end.eq.0 .and. nsup.gt.0) then
934 else if (nz_end.eq.0) then
937 if(me.eq.king.or..not.out1file)then
938 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
939 write (iout,*) 'IZ_SC=',iz_sc
941 c----------------------
944 if (.not.pdbref) then
945 call read_angles(inp,*38)
947 38 write (iout,'(a)') 'Error reading reference structure.'
949 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
950 stop 'Error reading reference structure'
954 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
964 call contact(.true.,ncont_ref,icont_ref,co)
966 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
968 if (constr_dist.gt.0) call read_dist_constr
969 write (iout,*) "After read_dist_constr nhpb",nhpb
970 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
972 if(me.eq.king.or..not.out1file)
973 & write (iout,*) 'Contact order:',co
975 if(me.eq.king.or..not.out1file)
976 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
979 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
981 if(me.eq.king.or..not.out1file)
982 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
983 & icont_ref(1,i),' ',
984 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
988 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
989 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
990 & modecalc.ne.10) then
991 C If input structure hasn't been supplied from the PDB file read or generate
993 if (iranconf.eq.0 .and. .not. extconf) then
994 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
995 & write (iout,'(a)') 'Initial geometry will be read in.'
997 read(inp,'(8f10.5)',end=36,err=36)
998 & ((c(l,k),l=1,3),k=1,nres),
999 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1000 write (iout,*) "Exit READ_CART"
1001 write (iout,'(8f10.5)')
1002 & ((c(l,k),l=1,3),k=1,nres),
1003 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1004 call int_from_cart1(.true.)
1005 write (iout,*) "Finish INT_TO_CART"
1008 dc(j,i)=c(j,i+1)-c(j,i)
1009 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1013 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1015 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1016 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1022 call read_angles(inp,*36)
1025 36 write (iout,'(a)') 'Error reading angle file.'
1027 call mpi_finalize( MPI_COMM_WORLD,IERR )
1029 stop 'Error reading angle file.'
1031 else if (extconf) then
1032 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1033 & write (iout,'(a)') 'Extended chain initial geometry.'
1035 theta(i)=90d0*deg2rad
1038 phi(i)=180d0*deg2rad
1041 alph(i)=110d0*deg2rad
1044 omeg(i)=-120d0*deg2rad
1045 if (itype(i).le.0) omeg(i)=-omeg(i)
1048 if(me.eq.king.or..not.out1file)
1049 & write (iout,'(a)') 'Random-generated initial geometry.'
1053 if (me.eq.king .or. fg_rank.eq.0 .and. (
1054 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1058 call gen_rand_conf(itmp,*30)
1060 30 write (iout,*) 'Failed to generate random conformation',
1061 & ', itrial=',itrial
1062 write (*,*) 'Processor:',me,
1063 & ' Failed to generate random conformation',
1073 write (iout,'(a,i3,a)') 'Processor:',me,
1074 & ' error in generating random conformation.'
1075 write (*,'(a,i3,a)') 'Processor:',me,
1076 & ' error in generating random conformation.'
1079 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1084 & ' error in generating random conformation.'
1089 elseif (modecalc.eq.4) then
1090 read (inp,'(a)') intinname
1091 open (intin,file=intinname,status='old',err=333)
1092 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1093 & write (iout,'(a)') 'intinname',intinname
1094 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1096 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1098 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1100 stop 'Error opening angle file.'
1104 C Generate distance constraints, if the PDB structure is to be regularized.
1105 if (nthread.gt.0) then
1106 call read_threadbase
1109 if (me.eq.king .or. .not. out1file)
1111 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1112 write (iout,'(/a,i3,a)')
1113 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1114 write (iout,'(20i4)') (iss(i),i=1,ns)
1116 write(iout,*)"Running with dynamic disulfide-bond formation"
1118 write (iout,'(/a/)') 'Pre-formed links are:'
1124 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1125 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1131 if (ns.gt.0.and.dyn_ss) then
1135 forcon(i-nss)=forcon(i)
1142 dyn_ss_mask(iss(i))=.true.
1145 if (i2ndstr.gt.0) call secstrp2dihc
1146 c call geom_to_var(nvar,x)
1147 c call etotal(energia(0))
1148 c call enerprint(energia(0))
1149 c call briefout(0,etot)
1151 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1152 cd write (iout,'(a)') 'Variable list:'
1153 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1155 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1156 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1157 & 'Processor',myrank,': end reading molecular data.'
1161 c--------------------------------------------------------------------------
1162 logical function seq_comp(itypea,itypeb,length)
1164 integer length,itypea(length),itypeb(length)
1167 if (itypea(i).ne.itypeb(i)) then
1175 c-----------------------------------------------------------------------------
1176 subroutine read_bridge
1177 C Read information about disulfide bridges.
1178 implicit real*8 (a-h,o-z)
1179 include 'DIMENSIONS'
1183 include 'COMMON.IOUNITS'
1184 include 'COMMON.GEO'
1185 include 'COMMON.VAR'
1186 include 'COMMON.INTERACT'
1187 include 'COMMON.LOCAL'
1188 include 'COMMON.NAMES'
1189 include 'COMMON.CHAIN'
1190 include 'COMMON.FFIELD'
1191 include 'COMMON.SBRIDGE'
1192 include 'COMMON.HEADER'
1193 include 'COMMON.CONTROL'
1194 include 'COMMON.DBASE'
1195 include 'COMMON.THREAD'
1196 include 'COMMON.TIME1'
1197 include 'COMMON.SETUP'
1198 C Read bridging residues.
1199 read (inp,*) ns,(iss(i),i=1,ns)
1201 if(me.eq.king.or..not.out1file)
1202 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1203 C Check whether the specified bridging residues are cystines.
1205 if (itype(iss(i)).ne.1) then
1206 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1207 & 'Do you REALLY think that the residue ',
1208 & restyp(itype(iss(i))),i,
1209 & ' can form a disulfide bridge?!!!'
1210 write (*,'(2a,i3,a)')
1211 & 'Do you REALLY think that the residue ',
1212 & restyp(itype(iss(i))),i,
1213 & ' can form a disulfide bridge?!!!'
1215 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1220 C Read preformed bridges.
1222 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1224 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1227 C Check if the residues involved in bridges are in the specified list of
1228 C bridging residues.
1231 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1232 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1233 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1234 & ' contains residues present in other pairs.'
1235 write (*,'(a,i3,a)') 'Disulfide pair',i,
1236 & ' contains residues present in other pairs.'
1238 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1244 if (ihpb(i).eq.iss(j)) goto 10
1246 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1249 if (jhpb(i).eq.iss(j)) goto 20
1251 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1257 ihpb(i)=ihpb(i)+nres
1258 jhpb(i)=jhpb(i)+nres
1264 c----------------------------------------------------------------------------
1265 subroutine read_x(kanal,*)
1266 implicit real*8 (a-h,o-z)
1267 include 'DIMENSIONS'
1268 include 'COMMON.GEO'
1269 include 'COMMON.VAR'
1270 include 'COMMON.CHAIN'
1271 include 'COMMON.IOUNITS'
1272 include 'COMMON.CONTROL'
1273 include 'COMMON.LOCAL'
1274 include 'COMMON.INTERACT'
1275 c Read coordinates from input
1277 read(kanal,'(8f10.5)',end=10,err=10)
1278 & ((c(l,k),l=1,3),k=1,nres),
1279 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1282 c(j,2*nres)=c(j,nres)
1284 call int_from_cart1(.false.)
1287 dc(j,i)=c(j,i+1)-c(j,i)
1288 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1292 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1294 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1295 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1303 c----------------------------------------------------------------------------
1304 subroutine read_threadbase
1305 implicit real*8 (a-h,o-z)
1306 include 'DIMENSIONS'
1307 include 'COMMON.IOUNITS'
1308 include 'COMMON.GEO'
1309 include 'COMMON.VAR'
1310 include 'COMMON.INTERACT'
1311 include 'COMMON.LOCAL'
1312 include 'COMMON.NAMES'
1313 include 'COMMON.CHAIN'
1314 include 'COMMON.FFIELD'
1315 include 'COMMON.SBRIDGE'
1316 include 'COMMON.HEADER'
1317 include 'COMMON.CONTROL'
1318 include 'COMMON.DBASE'
1319 include 'COMMON.THREAD'
1320 include 'COMMON.TIME1'
1321 C Read pattern database for threading.
1322 read (icbase,*) nseq
1324 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1325 & nres_base(2,i),nres_base(3,i)
1326 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1328 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1329 c & nres_base(2,i),nres_base(3,i)
1330 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1334 if (weidis.eq.0.0D0) weidis=0.1D0
1343 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1344 write (iout,'(a,i5)') 'nexcl: ',nexcl
1345 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1348 c------------------------------------------------------------------------------
1349 subroutine setup_var
1350 implicit real*8 (a-h,o-z)
1351 include 'DIMENSIONS'
1352 include 'COMMON.IOUNITS'
1353 include 'COMMON.GEO'
1354 include 'COMMON.VAR'
1355 include 'COMMON.INTERACT'
1356 include 'COMMON.LOCAL'
1357 include 'COMMON.NAMES'
1358 include 'COMMON.CHAIN'
1359 include 'COMMON.FFIELD'
1360 include 'COMMON.SBRIDGE'
1361 include 'COMMON.HEADER'
1362 include 'COMMON.CONTROL'
1363 include 'COMMON.DBASE'
1364 include 'COMMON.THREAD'
1365 include 'COMMON.TIME1'
1366 C Set up variable list.
1372 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1374 ialph(i,1)=nvar+nside
1378 if (indphi.gt.0) then
1380 else if (indback.gt.0) then
1385 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1388 c----------------------------------------------------------------------------
1389 subroutine gen_dist_constr
1390 C Generate CA distance constraints.
1391 implicit real*8 (a-h,o-z)
1392 include 'DIMENSIONS'
1393 include 'COMMON.IOUNITS'
1394 include 'COMMON.GEO'
1395 include 'COMMON.VAR'
1396 include 'COMMON.INTERACT'
1397 include 'COMMON.LOCAL'
1398 include 'COMMON.NAMES'
1399 include 'COMMON.CHAIN'
1400 include 'COMMON.FFIELD'
1401 include 'COMMON.SBRIDGE'
1402 include 'COMMON.HEADER'
1403 include 'COMMON.CONTROL'
1404 include 'COMMON.DBASE'
1405 include 'COMMON.THREAD'
1406 include 'COMMON.TIME1'
1407 dimension itype_pdb(maxres)
1408 common /pizda/ itype_pdb
1410 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1411 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1412 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1414 do i=nstart_sup,nstart_sup+nsup-1
1415 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1416 cd & ' seq_pdb', restyp(itype_pdb(i))
1417 do j=i+2,nstart_sup+nsup-1
1419 ihpb(nhpb)=i+nstart_seq-nstart_sup
1420 jhpb(nhpb)=j+nstart_seq-nstart_sup
1422 dhpb(nhpb)=dist(i,j)
1425 cd write (iout,'(a)') 'Distance constraints:'
1430 cd if (ii.gt.nres) then
1435 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1436 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1437 cd & dhpb(i),forcon(i)
1441 c----------------------------------------------------------------------------
1443 implicit real*8 (a-h,o-z)
1444 include 'DIMENSIONS'
1445 include 'COMMON.MAP'
1446 include 'COMMON.IOUNITS'
1447 character*3 angid(4) /'THE','PHI','ALP','OME'/
1448 character*80 mapcard,ucase
1450 read (inp,'(a)') mapcard
1451 mapcard=ucase(mapcard)
1452 if (index(mapcard,'PHI').gt.0) then
1454 else if (index(mapcard,'THE').gt.0) then
1456 else if (index(mapcard,'ALP').gt.0) then
1458 else if (index(mapcard,'OME').gt.0) then
1461 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1462 stop 'Error - illegal variable spec in MAP card.'
1464 call readi (mapcard,'RES1',res1(imap),0)
1465 call readi (mapcard,'RES2',res2(imap),0)
1466 if (res1(imap).eq.0) then
1467 res1(imap)=res2(imap)
1468 else if (res2(imap).eq.0) then
1469 res2(imap)=res1(imap)
1471 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1473 & 'Error - illegal definition of variable group in MAP.'
1474 stop 'Error - illegal definition of variable group in MAP.'
1476 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1477 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1478 call readi(mapcard,'NSTEP',nstep(imap),0)
1479 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1481 & 'Illegal boundary and/or step size specification in MAP.'
1482 stop 'Illegal boundary and/or step size specification in MAP.'
1487 c----------------------------------------------------------------------------
1489 implicit real*8 (a-h,o-z)
1490 include 'DIMENSIONS'
1491 include 'COMMON.IOUNITS'
1492 include 'COMMON.GEO'
1493 include 'COMMON.CSA'
1494 include 'COMMON.BANK'
1495 include 'COMMON.CONTROL'
1497 character*620 mcmcard
1498 call card_concat(mcmcard)
1500 call readi(mcmcard,'NCONF',nconf,50)
1501 call readi(mcmcard,'NADD',nadd,0)
1502 call readi(mcmcard,'JSTART',jstart,1)
1503 call readi(mcmcard,'JEND',jend,1)
1504 call readi(mcmcard,'NSTMAX',nstmax,500000)
1505 call readi(mcmcard,'N0',n0,1)
1506 call readi(mcmcard,'N1',n1,6)
1507 call readi(mcmcard,'N2',n2,4)
1508 call readi(mcmcard,'N3',n3,0)
1509 call readi(mcmcard,'N4',n4,0)
1510 call readi(mcmcard,'N5',n5,0)
1511 call readi(mcmcard,'N6',n6,10)
1512 call readi(mcmcard,'N7',n7,0)
1513 call readi(mcmcard,'N8',n8,0)
1514 call readi(mcmcard,'N9',n9,0)
1515 call readi(mcmcard,'N14',n14,0)
1516 call readi(mcmcard,'N15',n15,0)
1517 call readi(mcmcard,'N16',n16,0)
1518 call readi(mcmcard,'N17',n17,0)
1519 call readi(mcmcard,'N18',n18,0)
1521 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1523 call readi(mcmcard,'NDIFF',ndiff,2)
1524 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1525 call readi(mcmcard,'IS1',is1,1)
1526 call readi(mcmcard,'IS2',is2,8)
1527 call readi(mcmcard,'NRAN0',nran0,4)
1528 call readi(mcmcard,'NRAN1',nran1,2)
1529 call readi(mcmcard,'IRR',irr,1)
1530 call readi(mcmcard,'NSEED',nseed,20)
1531 call readi(mcmcard,'NTOTAL',ntotal,10000)
1532 call reada(mcmcard,'CUT1',cut1,2.0d0)
1533 call reada(mcmcard,'CUT2',cut2,5.0d0)
1534 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1535 call readi(mcmcard,'ICMAX',icmax,3)
1536 call readi(mcmcard,'IRESTART',irestart,0)
1537 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1540 call reada(mcmcard,'DELE',dele,20.0d0)
1541 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1542 call readi(mcmcard,'IREF',iref,0)
1543 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1544 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1545 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1546 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1547 write (iout,*) "NCONF_IN",nconf_in
1550 c----------------------------------------------------------------------------
1551 cfmc subroutine mcmfread
1552 cfmc implicit real*8 (a-h,o-z)
1553 cfmc include 'DIMENSIONS'
1554 cfmc include 'COMMON.MCMF'
1555 cfmc include 'COMMON.IOUNITS'
1556 cfmc include 'COMMON.GEO'
1557 cfmc character*80 ucase
1558 cfmc character*620 mcmcard
1559 cfmc call card_concat(mcmcard)
1561 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1562 cfmc write(iout,*)'MAXRANT=',maxrant
1563 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1564 cfmc write(iout,*)'MAXFAM=',maxfam
1565 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1566 cfmc write(iout,*)'NNET1=',nnet1
1567 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1568 cfmc write(iout,*)'NNET2=',nnet2
1569 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1570 cfmc write(iout,*)'NNET3=',nnet3
1571 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1572 cfmc write(iout,*)'ILASTT=',ilastt
1573 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1574 cfmc write(iout,*)'MAXSTR=',maxstr
1575 cfmc maxstr_f=maxstr/maxfam
1576 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1577 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1578 cfmc write(iout,*)'NMCMF=',nmcmf
1579 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1580 cfmc write(iout,*)'IFOCUS=',ifocus
1581 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1582 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1583 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1584 cfmc write(iout,*)'INTPRT=',intprt
1585 cfmc call readi(mcmcard,'IPRT',iprt,100)
1586 cfmc write(iout,*)'IPRT=',iprt
1587 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1588 cfmc write(iout,*)'IMAXTR=',imaxtr
1589 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1590 cfmc write(iout,*)'MAXEVEN=',maxeven
1591 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1592 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1593 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1594 cfmc write(iout,*)'INIMIN=',inimin
1595 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1596 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1597 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1598 cfmc write(iout,*)'NTHREAD=',nthread
1599 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1600 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1601 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1602 cfmc write(iout,*)'MAXPERT=',maxpert
1603 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1604 cfmc write(iout,*)'IRMSD=',irmsd
1605 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1606 cfmc write(iout,*)'DENEMIN=',denemin
1607 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1608 cfmc write(iout,*)'RCUT1S=',rcut1s
1609 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1610 cfmc write(iout,*)'RCUT1E=',rcut1e
1611 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1612 cfmc write(iout,*)'RCUT2S=',rcut2s
1613 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1614 cfmc write(iout,*)'RCUT2E=',rcut2e
1615 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1616 cfmc write(iout,*)'DPERT1=',d_pert1
1617 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1618 cfmc write(iout,*)'DPERT1A=',d_pert1a
1619 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1620 cfmc write(iout,*)'DPERT2=',d_pert2
1621 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1622 cfmc write(iout,*)'DPERT2A=',d_pert2a
1623 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1624 cfmc write(iout,*)'DPERT2B=',d_pert2b
1625 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1626 cfmc write(iout,*)'DPERT2C=',d_pert2c
1627 cfmc d_pert1=deg2rad*d_pert1
1628 cfmc d_pert1a=deg2rad*d_pert1a
1629 cfmc d_pert2=deg2rad*d_pert2
1630 cfmc d_pert2a=deg2rad*d_pert2a
1631 cfmc d_pert2b=deg2rad*d_pert2b
1632 cfmc d_pert2c=deg2rad*d_pert2c
1633 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1634 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1635 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1636 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1637 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1638 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1639 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1640 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1641 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1642 cfmc write(iout,*)'RCUTINI=',rcutini
1643 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1644 cfmc write(iout,*)'GRAT=',grat
1645 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1646 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1650 c----------------------------------------------------------------------------
1652 implicit real*8 (a-h,o-z)
1653 include 'DIMENSIONS'
1654 include 'COMMON.MCM'
1655 include 'COMMON.MCE'
1656 include 'COMMON.IOUNITS'
1658 character*320 mcmcard
1659 call card_concat(mcmcard)
1660 call readi(mcmcard,'MAXACC',maxacc,100)
1661 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1662 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1663 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1664 call readi(mcmcard,'MAXREPM',maxrepm,200)
1665 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1666 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1667 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1668 call reada(mcmcard,'E_UP',e_up,5.0D0)
1669 call reada(mcmcard,'DELTE',delte,0.1D0)
1670 call readi(mcmcard,'NSWEEP',nsweep,5)
1671 call readi(mcmcard,'NSTEPH',nsteph,0)
1672 call readi(mcmcard,'NSTEPC',nstepc,0)
1673 call reada(mcmcard,'TMIN',tmin,298.0D0)
1674 call reada(mcmcard,'TMAX',tmax,298.0D0)
1675 call readi(mcmcard,'NWINDOW',nwindow,0)
1676 call readi(mcmcard,'PRINT_MC',print_mc,0)
1677 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1678 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1679 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1680 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1681 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1682 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1683 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1684 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1685 if (nwindow.gt.0) then
1686 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1688 winlen(i)=winend(i)-winstart(i)+1
1691 if (tmax.lt.tmin) tmax=tmin
1692 if (tmax.eq.tmin) then
1696 if (nstepc.gt.0 .and. nsteph.gt.0) then
1697 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1698 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1700 C Probabilities of different move types
1701 sumpro_type(0)=0.0D0
1702 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1703 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1704 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1705 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1706 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1707 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1708 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1710 print *,'i',i,' sumprotype',sumpro_type(i)
1711 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1712 print *,'i',i,' sumprotype',sumpro_type(i)
1716 c----------------------------------------------------------------------------
1717 subroutine read_minim
1718 implicit real*8 (a-h,o-z)
1719 include 'DIMENSIONS'
1720 include 'COMMON.MINIM'
1721 include 'COMMON.IOUNITS'
1723 character*320 minimcard
1724 call card_concat(minimcard)
1725 call readi(minimcard,'MAXMIN',maxmin,2000)
1726 call readi(minimcard,'MAXFUN',maxfun,5000)
1727 call readi(minimcard,'MINMIN',minmin,maxmin)
1728 call readi(minimcard,'MINFUN',minfun,maxmin)
1729 call reada(minimcard,'TOLF',tolf,1.0D-2)
1730 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1731 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1732 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1733 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1734 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1735 & 'Options in energy minimization:'
1736 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1737 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1738 & 'MinMin:',MinMin,' MinFun:',MinFun,
1739 & ' TolF:',TolF,' RTolF:',RTolF
1742 c----------------------------------------------------------------------------
1743 subroutine read_angles(kanal,*)
1744 implicit real*8 (a-h,o-z)
1745 include 'DIMENSIONS'
1746 include 'COMMON.GEO'
1747 include 'COMMON.VAR'
1748 include 'COMMON.CHAIN'
1749 include 'COMMON.IOUNITS'
1750 include 'COMMON.CONTROL'
1751 c Read angles from input
1753 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1754 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1755 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1756 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1759 c 9/7/01 avoid 180 deg valence angle
1760 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1762 theta(i)=deg2rad*theta(i)
1763 phi(i)=deg2rad*phi(i)
1764 alph(i)=deg2rad*alph(i)
1765 omeg(i)=deg2rad*omeg(i)
1770 c----------------------------------------------------------------------------
1771 subroutine reada(rekord,lancuch,wartosc,default)
1773 character*(*) rekord,lancuch
1774 double precision wartosc,default
1777 iread=index(rekord,lancuch)
1778 if (iread.eq.0) then
1782 iread=iread+ilen(lancuch)+1
1783 read (rekord(iread:),*,err=10,end=10) wartosc
1788 c----------------------------------------------------------------------------
1789 subroutine readi(rekord,lancuch,wartosc,default)
1791 character*(*) rekord,lancuch
1792 integer wartosc,default
1795 iread=index(rekord,lancuch)
1796 if (iread.eq.0) then
1800 iread=iread+ilen(lancuch)+1
1801 read (rekord(iread:),*,err=10,end=10) wartosc
1806 c----------------------------------------------------------------------------
1807 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1810 integer tablica(dim),default
1811 character*(*) rekord,lancuch
1818 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1819 if (iread.eq.0) return
1820 iread=iread+ilen(lancuch)+1
1821 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1824 c----------------------------------------------------------------------------
1825 subroutine multreada(rekord,lancuch,tablica,dim,default)
1828 double precision tablica(dim),default
1829 character*(*) rekord,lancuch
1836 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1837 if (iread.eq.0) return
1838 iread=iread+ilen(lancuch)+1
1839 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1842 c----------------------------------------------------------------------------
1843 subroutine openunits
1844 implicit real*8 (a-h,o-z)
1845 include 'DIMENSIONS'
1848 character*16 form,nodename
1851 include 'COMMON.SETUP'
1852 include 'COMMON.IOUNITS'
1854 include 'COMMON.CONTROL'
1855 integer lenpre,lenpot,ilen,lentmp
1857 character*3 out1file_text,ucase
1860 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1861 call getenv_loc("PREFIX",prefix)
1863 call getenv_loc("POT",pot)
1864 call getenv_loc("DIRTMP",tmpdir)
1865 call getenv_loc("CURDIR",curdir)
1866 call getenv_loc("OUT1FILE",out1file_text)
1867 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1868 out1file_text=ucase(out1file_text)
1869 if (out1file_text(1:1).eq."Y") then
1872 out1file=fg_rank.gt.0
1877 if (lentmp.gt.0) then
1878 write (*,'(80(1h!))')
1879 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1880 write (*,'(80(1h!))')
1881 write (*,*)"All output files will be on node /tmp directory."
1883 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1884 if (me.eq.king) then
1885 write (*,*) "The master node is ",nodename
1886 else if (fg_rank.eq.0) then
1887 write (*,*) "I am the CG slave node ",nodename
1889 write (*,*) "I am the FG slave node ",nodename
1892 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1893 lenpre = lentmp+lenpre+1
1895 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1896 C Get the names and open the input files
1897 #if defined(WINIFL) || defined(WINPGI)
1898 open(1,file=pref_orig(:ilen(pref_orig))//
1899 & '.inp',status='old',readonly,shared)
1900 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1901 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1902 C Get parameter filenames and open the parameter files.
1903 call getenv_loc('BONDPAR',bondname)
1904 open (ibond,file=bondname,status='old',readonly,shared)
1905 call getenv_loc('THETPAR',thetname)
1906 open (ithep,file=thetname,status='old',readonly,shared)
1907 call getenv_loc('ROTPAR',rotname)
1908 open (irotam,file=rotname,status='old',readonly,shared)
1909 call getenv_loc('TORPAR',torname)
1910 open (itorp,file=torname,status='old',readonly,shared)
1911 call getenv_loc('TORDPAR',tordname)
1912 open (itordp,file=tordname,status='old',readonly,shared)
1913 call getenv_loc('FOURIER',fouriername)
1914 open (ifourier,file=fouriername,status='old',readonly,shared)
1915 call getenv_loc('ELEPAR',elename)
1916 open (ielep,file=elename,status='old',readonly,shared)
1917 call getenv_loc('SIDEPAR',sidename)
1918 open (isidep,file=sidename,status='old',readonly,shared)
1919 #elif (defined CRAY) || (defined AIX)
1920 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1922 c print *,"Processor",myrank," opened file 1"
1923 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1924 c print *,"Processor",myrank," opened file 9"
1925 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1926 C Get parameter filenames and open the parameter files.
1927 call getenv_loc('BONDPAR',bondname)
1928 open (ibond,file=bondname,status='old',action='read')
1929 c print *,"Processor",myrank," opened file IBOND"
1930 call getenv_loc('THETPAR',thetname)
1931 open (ithep,file=thetname,status='old',action='read')
1932 c print *,"Processor",myrank," opened file ITHEP"
1933 call getenv_loc('ROTPAR',rotname)
1934 open (irotam,file=rotname,status='old',action='read')
1935 c print *,"Processor",myrank," opened file IROTAM"
1936 call getenv_loc('TORPAR',torname)
1937 open (itorp,file=torname,status='old',action='read')
1938 c print *,"Processor",myrank," opened file ITORP"
1939 call getenv_loc('TORDPAR',tordname)
1940 open (itordp,file=tordname,status='old',action='read')
1941 c print *,"Processor",myrank," opened file ITORDP"
1942 call getenv_loc('SCCORPAR',sccorname)
1943 open (isccor,file=sccorname,status='old',action='read')
1944 c print *,"Processor",myrank," opened file ISCCOR"
1945 call getenv_loc('FOURIER',fouriername)
1946 open (ifourier,file=fouriername,status='old',action='read')
1947 c print *,"Processor",myrank," opened file IFOURIER"
1948 call getenv_loc('ELEPAR',elename)
1949 open (ielep,file=elename,status='old',action='read')
1950 c print *,"Processor",myrank," opened file IELEP"
1951 call getenv_loc('SIDEPAR',sidename)
1952 open (isidep,file=sidename,status='old',action='read')
1953 c print *,"Processor",myrank," opened file ISIDEP"
1954 c print *,"Processor",myrank," opened parameter files"
1956 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1957 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1958 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1959 C Get parameter filenames and open the parameter files.
1960 call getenv_loc('BONDPAR',bondname)
1961 open (ibond,file=bondname,status='old')
1962 call getenv_loc('THETPAR',thetname)
1963 open (ithep,file=thetname,status='old')
1964 call getenv_loc('ROTPAR',rotname)
1965 open (irotam,file=rotname,status='old')
1966 call getenv_loc('TORPAR',torname)
1967 open (itorp,file=torname,status='old')
1968 call getenv_loc('TORDPAR',tordname)
1969 open (itordp,file=tordname,status='old')
1970 call getenv_loc('SCCORPAR',sccorname)
1971 open (isccor,file=sccorname,status='old')
1972 call getenv_loc('FOURIER',fouriername)
1973 open (ifourier,file=fouriername,status='old')
1974 call getenv_loc('ELEPAR',elename)
1975 open (ielep,file=elename,status='old')
1976 call getenv_loc('SIDEPAR',sidename)
1977 open (isidep,file=sidename,status='old')
1979 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1981 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1982 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1983 C Get parameter filenames and open the parameter files.
1984 call getenv_loc('BONDPAR',bondname)
1985 open (ibond,file=bondname,status='old',readonly)
1986 call getenv_loc('THETPAR',thetname)
1987 open (ithep,file=thetname,status='old',readonly)
1988 call getenv_loc('ROTPAR',rotname)
1989 open (irotam,file=rotname,status='old',readonly)
1990 call getenv_loc('TORPAR',torname)
1991 open (itorp,file=torname,status='old',readonly)
1992 call getenv_loc('TORDPAR',tordname)
1993 open (itordp,file=tordname,status='old',readonly)
1994 call getenv_loc('SCCORPAR',sccorname)
1995 open (isccor,file=sccorname,status='old',readonly)
1997 call getenv_loc('THETPARPDB',thetname_pdb)
1998 print *,"thetname_pdb ",thetname_pdb
1999 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2000 print *,ithep_pdb," opened"
2002 call getenv_loc('FOURIER',fouriername)
2003 open (ifourier,file=fouriername,status='old',readonly)
2004 call getenv_loc('ELEPAR',elename)
2005 open (ielep,file=elename,status='old',readonly)
2006 call getenv_loc('SIDEPAR',sidename)
2007 open (isidep,file=sidename,status='old',readonly)
2008 call getenv_loc('LIPTRANPAR',liptranname)
2009 open (iliptranpar,file=liptranname,status='old',action='read')
2011 call getenv_loc('ROTPARPDB',rotname_pdb)
2012 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2017 C 8/9/01 In the newest version SCp interaction constants are read from a file
2018 C Use -DOLDSCP to use hard-coded constants instead.
2020 call getenv_loc('SCPPAR',scpname)
2021 #if defined(WINIFL) || defined(WINPGI)
2022 open (iscpp,file=scpname,status='old',readonly,shared)
2023 #elif (defined CRAY) || (defined AIX)
2024 open (iscpp,file=scpname,status='old',action='read')
2026 open (iscpp,file=scpname,status='old')
2028 open (iscpp,file=scpname,status='old',readonly)
2031 call getenv_loc('PATTERN',patname)
2032 #if defined(WINIFL) || defined(WINPGI)
2033 open (icbase,file=patname,status='old',readonly,shared)
2034 #elif (defined CRAY) || (defined AIX)
2035 open (icbase,file=patname,status='old',action='read')
2037 open (icbase,file=patname,status='old')
2039 open (icbase,file=patname,status='old',readonly)
2042 C Open output file only for CG processes
2043 c print *,"Processor",myrank," fg_rank",fg_rank
2044 if (fg_rank.eq.0) then
2046 if (nodes.eq.1) then
2049 npos = dlog10(dfloat(nodes-1))+1
2051 if (npos.lt.3) npos=3
2052 write (liczba,'(i1)') npos
2053 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2055 write (liczba,form) me
2056 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2057 & liczba(:ilen(liczba))
2058 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2060 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2062 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2063 & liczba(:ilen(liczba))//'.mol2'
2064 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2065 & liczba(:ilen(liczba))//'.stat'
2067 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2068 & //liczba(:ilen(liczba))//'.stat')
2069 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2072 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2073 & liczba(:ilen(liczba))//'.const'
2078 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2079 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2080 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2081 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2082 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2084 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2086 rest2name=prefix(:ilen(prefix))//'.rst'
2088 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2091 #if defined(AIX) || defined(PGI)
2092 if (me.eq.king .or. .not. out1file)
2093 & open(iout,file=outname,status='unknown')
2095 if (fg_rank.gt.0) then
2096 write (liczba,'(i3.3)') myrank/nfgtasks
2097 write (ll,'(bz,i3.3)') fg_rank
2098 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2103 open(igeom,file=intname,status='unknown',position='append')
2104 open(ipdb,file=pdbname,status='unknown')
2105 open(imol2,file=mol2name,status='unknown')
2106 open(istat,file=statname,status='unknown',position='append')
2108 c1out open(iout,file=outname,status='unknown')
2111 if (me.eq.king .or. .not.out1file)
2112 & open(iout,file=outname,status='unknown')
2114 if (fg_rank.gt.0) then
2115 write (liczba,'(i3.3)') myrank/nfgtasks
2116 write (ll,'(bz,i3.3)') fg_rank
2117 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2122 open(igeom,file=intname,status='unknown',access='append')
2123 open(ipdb,file=pdbname,status='unknown')
2124 open(imol2,file=mol2name,status='unknown')
2125 open(istat,file=statname,status='unknown',access='append')
2127 c1out open(iout,file=outname,status='unknown')
2130 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2131 csa_seed=prefix(:lenpre)//'.CSA.seed'
2132 csa_history=prefix(:lenpre)//'.CSA.history'
2133 csa_bank=prefix(:lenpre)//'.CSA.bank'
2134 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2135 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2136 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2137 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2138 csa_int=prefix(:lenpre)//'.int'
2139 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2140 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2141 csa_in=prefix(:lenpre)//'.CSA.in'
2142 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2145 write (iout,'(80(1h-))')
2146 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2147 write (iout,'(80(1h-))')
2148 write (iout,*) "Input file : ",
2149 & pref_orig(:ilen(pref_orig))//'.inp'
2150 write (iout,*) "Output file : ",
2151 & outname(:ilen(outname))
2153 write (iout,*) "Sidechain potential file : ",
2154 & sidename(:ilen(sidename))
2156 write (iout,*) "SCp potential file : ",
2157 & scpname(:ilen(scpname))
2159 write (iout,*) "Electrostatic potential file : ",
2160 & elename(:ilen(elename))
2161 write (iout,*) "Cumulant coefficient file : ",
2162 & fouriername(:ilen(fouriername))
2163 write (iout,*) "Torsional parameter file : ",
2164 & torname(:ilen(torname))
2165 write (iout,*) "Double torsional parameter file : ",
2166 & tordname(:ilen(tordname))
2167 write (iout,*) "SCCOR parameter file : ",
2168 & sccorname(:ilen(sccorname))
2169 write (iout,*) "Bond & inertia constant file : ",
2170 & bondname(:ilen(bondname))
2171 write (iout,*) "Bending parameter file : ",
2172 & thetname(:ilen(thetname))
2173 write (iout,*) "Rotamer parameter file : ",
2174 & rotname(:ilen(rotname))
2175 write (iout,*) "Thetpdb parameter file : ",
2176 & thetname_pdb(:ilen(thetname_pdb))
2177 write (iout,*) "Threading database : ",
2178 & patname(:ilen(patname))
2180 &write (iout,*)" DIRTMP : ",
2182 write (iout,'(80(1h-))')
2186 c----------------------------------------------------------------------------
2187 subroutine card_concat(card)
2188 implicit real*8 (a-h,o-z)
2189 include 'DIMENSIONS'
2190 include 'COMMON.IOUNITS'
2192 character*80 karta,ucase
2194 read (inp,'(a)') karta
2197 do while (karta(80:80).eq.'&')
2198 card=card(:ilen(card)+1)//karta(:79)
2199 read (inp,'(a)') karta
2202 card=card(:ilen(card)+1)//karta
2205 c----------------------------------------------------------------------------------
2207 implicit real*8 (a-h,o-z)
2208 include 'DIMENSIONS'
2209 include 'COMMON.CHAIN'
2210 include 'COMMON.IOUNITS'
2212 open(irest2,file=rest2name,status='unknown')
2213 read(irest2,*) totT,EK,potE,totE,t_bath
2216 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2219 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2222 read (irest2,*) iset
2227 c---------------------------------------------------------------------------------
2228 subroutine read_fragments
2229 implicit real*8 (a-h,o-z)
2230 include 'DIMENSIONS'
2234 include 'COMMON.SETUP'
2235 include 'COMMON.CHAIN'
2236 include 'COMMON.IOUNITS'
2238 include 'COMMON.CONTROL'
2239 read(inp,*) nset,nfrag,npair,nfrag_back
2240 if(me.eq.king.or..not.out1file)
2241 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2242 & " nfrag_back",nfrag_back
2244 read(inp,*) mset(iset)
2246 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2248 if(me.eq.king.or..not.out1file)
2249 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2250 & ifrag(2,i,iset), qinfrag(i,iset)
2253 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2255 if(me.eq.king.or..not.out1file)
2256 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2257 & ipair(2,i,iset), qinpair(i,iset)
2260 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2261 & wfrag_back(3,i,iset),
2262 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2263 if(me.eq.king.or..not.out1file)
2264 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2265 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2270 C---------------------------------------------------------------------------
2271 subroutine read_afminp
2272 implicit real*8 (a-h,o-z)
2273 include 'DIMENSIONS'
2277 include 'COMMON.SETUP'
2278 include 'COMMON.CONTROL'
2279 include 'COMMON.CHAIN'
2280 include 'COMMON.IOUNITS'
2281 include 'COMMON.SBRIDGE'
2282 character*320 afmcard
2284 call card_concat(afmcard)
2285 call readi(afmcard,"BEG",afmbeg,0)
2286 call readi(afmcard,"END",afmend,0)
2287 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2288 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2289 print *,'FORCE=' ,forceAFMconst
2290 CCCC NOW PROPERTIES FOR AFM
2293 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2295 distafminit=dsqrt(distafminit)
2296 print *,'initdist',distafminit
2300 c-------------------------------------------------------------------------------
2301 subroutine read_dist_constr
2302 implicit real*8 (a-h,o-z)
2303 include 'DIMENSIONS'
2307 include 'COMMON.SETUP'
2308 include 'COMMON.CONTROL'
2309 include 'COMMON.CHAIN'
2310 include 'COMMON.IOUNITS'
2311 include 'COMMON.SBRIDGE'
2312 integer ifrag_(2,100),ipair_(2,100)
2313 double precision wfrag_(100),wpair_(100)
2314 character*500 controlcard
2315 c write (iout,*) "Calling read_dist_constr"
2316 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2318 call card_concat(controlcard)
2319 call readi(controlcard,"NFRAG",nfrag_,0)
2320 call readi(controlcard,"NPAIR",npair_,0)
2321 call readi(controlcard,"NDIST",ndist_,0)
2322 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2323 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2324 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2325 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2326 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2327 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2328 c write (iout,*) "IFRAG"
2330 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2332 c write (iout,*) "IPAIR"
2334 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2338 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2339 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2340 & ifrag_(2,i)=nstart_sup+nsup-1
2341 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2343 if (wfrag_(i).gt.0.0d0) then
2344 do j=ifrag_(1,i),ifrag_(2,i)-1
2345 do k=j+1,ifrag_(2,i)
2346 c write (iout,*) "j",j," k",k
2348 if (constr_dist.eq.1) then
2353 forcon(nhpb)=wfrag_(i)
2354 else if (constr_dist.eq.2) then
2355 if (ddjk.le.dist_cut) then
2360 forcon(nhpb)=wfrag_(i)
2367 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2370 if (.not.out1file .or. me.eq.king)
2371 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2372 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2374 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2375 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2382 if (wpair_(i).gt.0.0d0) then
2390 do j=ifrag_(1,ii),ifrag_(2,ii)
2391 do k=ifrag_(1,jj),ifrag_(2,jj)
2395 forcon(nhpb)=wpair_(i)
2396 dhpb(nhpb)=dist(j,k)
2398 if (.not.out1file .or. me.eq.king)
2399 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2400 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2402 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2403 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2410 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2411 if (forcon(nhpb+1).gt.0.0d0) then
2413 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2415 if (.not.out1file .or. me.eq.king)
2416 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2417 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2419 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2420 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2427 c-------------------------------------------------------------------------------
2429 subroutine flush(iu)
2434 subroutine flush(iu)
2439 c------------------------------------------------------------------------------
2440 subroutine copy_to_tmp(source)
2441 include "DIMENSIONS"
2442 include "COMMON.IOUNITS"
2443 character*(*) source
2444 character* 256 tmpfile
2448 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2449 inquire(file=tmpfile,exist=ex)
2451 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2452 & " to temporary directory..."
2453 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2454 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2458 c------------------------------------------------------------------------------
2459 subroutine move_from_tmp(source)
2460 include "DIMENSIONS"
2461 include "COMMON.IOUNITS"
2462 character*(*) source
2465 write (*,*) "Moving ",source(:ilen(source)),
2466 & " from temporary directory to working directory"
2467 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2468 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2471 c------------------------------------------------------------------------------
2472 subroutine random_init(seed)
2474 C Initialize random number generator
2476 implicit real*8 (a-h,o-z)
2477 include 'DIMENSIONS'
2480 logical OKRandom, prng_restart
2482 integer iseed_array(4)
2484 include 'COMMON.IOUNITS'
2485 include 'COMMON.TIME1'
2486 include 'COMMON.THREAD'
2487 include 'COMMON.SBRIDGE'
2488 include 'COMMON.CONTROL'
2489 include 'COMMON.MCM'
2490 include 'COMMON.MAP'
2491 include 'COMMON.HEADER'
2492 include 'COMMON.CSA'
2493 include 'COMMON.CHAIN'
2494 include 'COMMON.MUCA'
2496 include 'COMMON.FFIELD'
2497 include 'COMMON.SETUP'
2498 iseed=-dint(dabs(seed))
2499 if (iseed.eq.0) then
2500 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2501 & 'Random seed undefined. The program will stop.'
2502 write (*,'(/80(1h*)/20x,a/80(1h*))')
2503 & 'Random seed undefined. The program will stop.'
2505 call mpi_finalize(mpi_comm_world,ierr)
2507 stop 'Bad random seed.'
2510 if (fg_rank.eq.0) then
2514 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2515 OKRandom = prng_restart(me,iseed)
2518 tmp=65536.0d0**(4-i)
2519 iseed_array(i) = dint(seed/tmp)
2520 seed=seed-iseed_array(i)*tmp
2523 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2524 & (iseed_array(i),i=1,4)
2525 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2526 & (iseed_array(i),i=1,4)
2527 OKRandom = prng_restart(me,iseed_array)
2530 c r1 = prng_next(me)
2531 r1=ran_number(0.0D0,1.0D0)
2533 & write (iout,*) 'ran_num',r1
2534 if (r1.lt.0.0d0) OKRandom=.false.
2536 if (.not.OKRandom) then
2537 write (iout,*) 'PRNG IS NOT WORKING!!!'
2538 print *,'PRNG IS NOT WORKING!!!'
2541 call mpi_abort(mpi_comm_world,error_msg,ierr)
2544 write (iout,*) 'too many processors for parallel prng'
2545 write (*,*) 'too many processors for parallel prng'
2553 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)