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 include 'COMMON.SHIELD'
85 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
86 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
88 character*320 controlcard
93 read (INP,'(a)') titel
94 call card_concat(controlcard)
95 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
96 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
97 call reada(controlcard,'SEED',seed,0.0D0)
98 call random_init(seed)
99 C Set up the time limit (caution! The time must be input in minutes!)
100 read_cart=index(controlcard,'READ_CART').gt.0
101 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
102 C this variable with_theta_constr is the variable which allow to read and execute the
103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
104 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
105 write (iout,*) "with_theta_constr ",with_theta_constr
106 call readi(controlcard,'SYM',symetr,1)
107 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
108 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
109 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
110 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
111 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
112 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
113 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
114 call reada(controlcard,'DRMS',drms,0.1D0)
115 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
117 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
118 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
119 write (iout,'(a,f10.1)')'DRMS = ',drms
120 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
121 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
123 call readi(controlcard,'NZ_START',nz_start,0)
124 call readi(controlcard,'NZ_END',nz_end,0)
125 call readi(controlcard,'IZ_SC',iz_sc,0)
127 safety = 60.0d0*safety
130 call reada(controlcard,"T_BATH",t_bath,300.0d0)
131 minim=(index(controlcard,'MINIMIZE').gt.0)
132 dccart=(index(controlcard,'CART').gt.0)
133 overlapsc=(index(controlcard,'OVERLAP').gt.0)
134 overlapsc=.not.overlapsc
135 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
136 searchsc=.not.searchsc
137 sideadd=(index(controlcard,'SIDEADD').gt.0)
138 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
139 outpdb=(index(controlcard,'PDBOUT').gt.0)
140 outmol2=(index(controlcard,'MOL2OUT').gt.0)
141 pdbref=(index(controlcard,'PDBREF').gt.0)
142 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
143 indpdb=index(controlcard,'PDBSTART')
144 extconf=(index(controlcard,'EXTCONF').gt.0)
145 AFMlog=(index(controlcard,'AFM'))
146 selfguide=(index(controlcard,'SELFGUIDE'))
147 print *,'AFMlog',AFMlog,selfguide,"KUPA"
148 call readi(controlcard,'TUBEMOD',tubelog,0)
149 write (iout,*) TUBElog,"TUBEMODE"
150 call readi(controlcard,'GENCONSTR',genconstr,0)
151 C write (iout,*) TUBElog,"TUBEMODE"
152 call readi(controlcard,'IPRINT',iprint,0)
153 C SHIELD keyword sets if the shielding effect of side-chains is used
154 C 0 denots no shielding is used all peptide are equally despite the
155 C solvent accesible area
156 C 1 the newly introduced function
157 C 2 reseved for further possible developement
158 call readi(controlcard,'SHIELD',shield_mode,0)
159 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
160 write(iout,*) "shield_mode",shield_mode
162 call readi(controlcard,'TORMODE',tor_mode,0)
163 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
164 write(iout,*) "torsional and valence angle mode",tor_mode
165 call readi(controlcard,'MAXGEN',maxgen,10000)
166 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
167 call readi(controlcard,"KDIAG",kdiag,0)
168 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
169 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
170 & write (iout,*) "RESCALE_MODE",rescale_mode
171 split_ene=index(controlcard,'SPLIT_ENE').gt.0
172 if (index(controlcard,'REGULAR').gt.0.0D0) then
173 call reada(controlcard,'WEIDIS',weidis,0.1D0)
177 if (index(controlcard,'CHECKGRAD').gt.0) then
179 if (index(controlcard,'CART').gt.0) then
181 elseif (index(controlcard,'CARINT').gt.0) then
186 call reada(controlcard,'DELTA',aincr,1.0d-7)
187 elseif (index(controlcard,'THREAD').gt.0) then
189 call readi(controlcard,'THREAD',nthread,0)
190 if (nthread.gt.0) then
191 call reada(controlcard,'WEIDIS',weidis,0.1D0)
194 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
195 stop 'Error termination in Read_Control.'
197 else if (index(controlcard,'MCMA').gt.0) then
199 else if (index(controlcard,'MCEE').gt.0) then
201 else if (index(controlcard,'MULTCONF').gt.0) then
203 else if (index(controlcard,'MAP').gt.0) then
205 call readi(controlcard,'MAP',nmap,0)
206 else if (index(controlcard,'CSA').gt.0) then
208 crc else if (index(controlcard,'ZSCORE').gt.0) then
210 crc ZSCORE is rm from UNRES, modecalc=9 is available
213 cfcm else if (index(controlcard,'MCMF').gt.0) then
215 else if (index(controlcard,'SOFTREG').gt.0) then
217 else if (index(controlcard,'CHECK_BOND').gt.0) then
219 else if (index(controlcard,'TEST').gt.0) then
221 else if (index(controlcard,'MD').gt.0) then
223 else if (index(controlcard,'RE ').gt.0) then
227 lmuca=index(controlcard,'MUCA').gt.0
228 call readi(controlcard,'MUCADYN',mucadyn,0)
229 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
230 if (lmuca .and. (me.eq.king .or. .not.out1file ))
232 write (iout,*) 'MUCADYN=',mucadyn
233 write (iout,*) 'MUCASMOOTH=',muca_smooth
236 iscode=index(controlcard,'ONE_LETTER')
237 indphi=index(controlcard,'PHI')
238 indback=index(controlcard,'BACK')
239 iranconf=index(controlcard,'RAND_CONF')
240 i2ndstr=index(controlcard,'USE_SEC_PRED')
241 gradout=index(controlcard,'GRADOUT').gt.0
242 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
243 C DISTCHAINMAX become obsolete for periodic boundry condition
244 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
245 C Reading the dimensions of box in x,y,z coordinates
246 call reada(controlcard,'BOXX',boxxsize,100.0d0)
247 call reada(controlcard,'BOXY',boxysize,100.0d0)
248 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
249 c Cutoff range for interactions
250 call reada(controlcard,"R_CUT",r_cut,15.0d0)
251 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
252 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
253 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
254 if (lipthick.gt.0.0d0) then
255 bordliptop=(boxzsize+lipthick)/2.0
256 bordlipbot=bordliptop-lipthick
258 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
259 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
260 buflipbot=bordlipbot+lipbufthick
261 bufliptop=bordliptop-lipbufthick
262 if ((lipbufthick*2.0d0).gt.lipthick)
263 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
265 write(iout,*) "bordliptop=",bordliptop
266 write(iout,*) "bordlipbot=",bordlipbot
267 write(iout,*) "bufliptop=",bufliptop
268 write(iout,*) "buflipbot=",buflipbot
269 write (iout,*) "SHIELD MODE",shield_mode
270 if (TUBElog.gt.0) then
271 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
272 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
273 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
274 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
275 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
276 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
277 buftubebot=bordtubebot+tubebufthick
278 buftubetop=bordtubetop-tubebufthick
280 if (shield_mode.gt.0) then
282 C VSolvSphere the volume of solving sphere
284 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
285 C there will be no distinction between proline peptide group and normal peptide
286 C group in case of shielding parameters
287 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
288 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
289 write (iout,*) VSolvSphere,VSolvSphere_div
290 C long axis of side chain
292 long_r_sidechain(i)=vbldsc0(1,i)
293 short_r_sidechain(i)=sigma0(i)
297 if (me.eq.king .or. .not.out1file )
298 & write (iout,*) "DISTCHAINMAX",distchainmax
300 if(me.eq.king.or..not.out1file)
301 & write (iout,'(2a)') diagmeth(kdiag),
302 & ' routine used to diagonalize matrices.'
305 c--------------------------------------------------------------------------
306 subroutine read_REMDpar
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.REMD'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SETUP'
327 character*320 controlcard
328 character*3200 controlcard1
329 integer iremd_m_total
331 if(me.eq.king.or..not.out1file)
332 & write (iout,*) "REMD setup"
334 call card_concat(controlcard)
335 call readi(controlcard,"NREP",nrep,3)
336 call readi(controlcard,"NSTEX",nstex,1000)
337 call reada(controlcard,"RETMIN",retmin,10.0d0)
338 call reada(controlcard,"RETMAX",retmax,1000.0d0)
339 mremdsync=(index(controlcard,'SYNC').gt.0)
340 call readi(controlcard,"NSYN",i_sync_step,100)
341 restart1file=(index(controlcard,'REST1FILE').gt.0)
342 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
343 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
344 if(max_cache_traj_use.gt.max_cache_traj)
345 & max_cache_traj_use=max_cache_traj
346 if(me.eq.king.or..not.out1file) then
347 cd if (traj1file) then
348 crc caching is in testing - NTWX is not ignored
349 cd write (iout,*) "NTWX value is ignored"
350 cd write (iout,*) " trajectory is stored to one file by master"
351 cd write (iout,*) " before exchange at NSTEX intervals"
353 write (iout,*) "NREP= ",nrep
354 write (iout,*) "NSTEX= ",nstex
355 write (iout,*) "SYNC= ",mremdsync
356 write (iout,*) "NSYN= ",i_sync_step
357 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
360 if (index(controlcard,'TLIST').gt.0) then
362 call card_concat(controlcard1)
363 read(controlcard1,*) (remd_t(i),i=1,nrep)
364 if(me.eq.king.or..not.out1file)
365 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
368 if (index(controlcard,'MLIST').gt.0) then
370 call card_concat(controlcard1)
371 read(controlcard1,*) (remd_m(i),i=1,nrep)
372 if(me.eq.king.or..not.out1file) then
373 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
376 iremd_m_total=iremd_m_total+remd_m(i)
378 write (iout,*) 'Total number of replicas ',iremd_m_total
381 if(me.eq.king.or..not.out1file)
382 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
385 c--------------------------------------------------------------------------
386 subroutine read_MDpar
390 implicit real*8 (a-h,o-z)
392 include 'COMMON.IOUNITS'
393 include 'COMMON.TIME1'
396 include 'COMMON.LANGEVIN'
398 include 'COMMON.LANGEVIN.lang0'
400 include 'COMMON.INTERACT'
401 include 'COMMON.NAMES'
403 include 'COMMON.SETUP'
404 include 'COMMON.CONTROL'
405 include 'COMMON.SPLITELE'
407 character*320 controlcard
409 call card_concat(controlcard)
410 call readi(controlcard,"NSTEP",n_timestep,1000000)
411 call readi(controlcard,"NTWE",ntwe,100)
412 call readi(controlcard,"NTWX",ntwx,1000)
413 call reada(controlcard,"DT",d_time,1.0d-1)
414 call reada(controlcard,"DVMAX",dvmax,2.0d1)
415 call reada(controlcard,"DAMAX",damax,1.0d1)
416 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
417 call readi(controlcard,"LANG",lang,0)
418 RESPA = index(controlcard,"RESPA") .gt. 0
419 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
420 ntime_split0=ntime_split
421 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
422 ntime_split0=ntime_split
423 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
424 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
425 rest = index(controlcard,"REST").gt.0
426 tbf = index(controlcard,"TBF").gt.0
427 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
428 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
429 tnh = index(controlcard,"NOSEHOOVER96").gt.0
430 if (RESPA.and.tnh)then
431 xiresp = index(controlcard,"XIRESP").gt.0
433 write(iout,*) "tnh",tnh
434 call reada(controlcard,"Q_NP",Q_np,0.1d0)
435 usampl = index(controlcard,"USAMPL").gt.0
437 mdpdb = index(controlcard,"MDPDB").gt.0
438 call reada(controlcard,"T_BATH",t_bath,300.0d0)
439 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
440 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
441 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
442 if (count_reset_moment.eq.0) count_reset_moment=1000000000
443 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
444 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
445 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
446 if (count_reset_vel.eq.0) count_reset_vel=1000000000
447 large = index(controlcard,"LARGE").gt.0
448 print_compon = index(controlcard,"PRINT_COMPON").gt.0
449 rattle = index(controlcard,"RATTLE").gt.0
450 c if performing umbrella sampling, fragments constrained are read from the fragment file
456 if(me.eq.king.or..not.out1file) then
458 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
460 write (iout,'(a)') "The units are:"
461 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
462 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
463 & " acceleration: angstrom/(48.9 fs)**2"
464 write (iout,'(a)') "energy: kcal/mol, temperature: K"
466 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
467 write (iout,'(a60,f10.5,a)')
468 & "Initial time step of numerical integration:",d_time,
470 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
472 write (iout,'(2a,i4,a)')
473 & "A-MTS algorithm used; initial time step for fast-varying",
474 & " short-range forces split into",ntime_split," steps."
475 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
476 & r_cut," lambda",rlamb
478 write (iout,'(2a,f10.5)')
479 & "Maximum acceleration threshold to reduce the time step",
480 & "/increase split number:",damax
481 write (iout,'(2a,f10.5)')
482 & "Maximum predicted energy drift to reduce the timestep",
483 & "/increase split number:",edriftmax
484 write (iout,'(a60,f10.5)')
485 & "Maximum velocity threshold to reduce velocities:",dvmax
486 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
487 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
488 if (rattle) write (iout,'(a60)')
489 & "Rattle algorithm used to constrain the virtual bonds"
493 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
494 call reada(controlcard,"RWAT",rwat,1.4d0)
495 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
496 surfarea=index(controlcard,"SURFAREA").gt.0
497 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
498 if(me.eq.king.or..not.out1file)then
499 write (iout,'(/a,$)') "Langevin dynamics calculation"
502 & " with direct integration of Langevin equations"
503 else if (lang.eq.2) then
504 write (iout,'(a/)') " with TINKER stochasic MD integrator"
505 else if (lang.eq.3) then
506 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
507 else if (lang.eq.4) then
508 write (iout,'(a/)') " in overdamped mode"
510 write (iout,'(//a,i5)')
511 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
514 write (iout,'(a60,f10.5)') "Temperature:",t_bath
515 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
516 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
517 write (iout,'(a60,f10.5)')
518 & "Scaling factor of the friction forces:",scal_fric
519 if (surfarea) write (iout,'(2a,i10,a)')
520 & "Friction coefficients will be scaled by solvent-accessible",
521 & " surface area every",reset_fricmat," steps."
523 c Calculate friction coefficients and bounds of stochastic forces
524 eta=6*pi*cPoise*etawat
525 if(me.eq.king.or..not.out1file)
526 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
528 gamp=scal_fric*(pstok+rwat)*eta
529 stdfp=dsqrt(2*Rb*t_bath/d_time)
531 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
532 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
534 if(me.eq.king.or..not.out1file)then
535 write (iout,'(/2a/)')
536 & "Radii of site types and friction coefficients and std's of",
537 & " stochastic forces of fully exposed sites"
538 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
540 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
541 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
545 if(me.eq.king.or..not.out1file)then
546 write (iout,'(a)') "Berendsen bath calculation"
547 write (iout,'(a60,f10.5)') "Temperature:",t_bath
548 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
550 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
551 & count_reset_moment," steps"
553 & write (iout,'(a,i10,a)')
554 & "Velocities will be reset at random every",count_reset_vel,
557 else if (tnp .or. tnp1 .or. tnh) then
558 if (tnp .or. tnp1) then
559 write (iout,'(a)') "Nose-Poincare bath calculation"
560 if (tnp) write (iout,'(a)')
561 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
562 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
564 write (iout,'(a)') "Nose-Hoover bath calculation"
565 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
575 WDTI(i) = 1.0*d_time/nresn
582 write (iout,'(a)') "NVT-XI-RESPA algorithm"
584 write (iout,'(a)') "NVT-XO-RESPA algorithm"
587 WDTIi(i) = 1.0*d_time/nresn/ntime_split
595 if(me.eq.king.or..not.out1file)
596 & write (iout,'(a31)') "Microcanonical mode calculation"
598 if(me.eq.king.or..not.out1file)then
599 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
601 write(iout,*) "MD running with constraints."
602 write(iout,*) "Equilibration time ", eq_time, " mtus."
603 write(iout,*) "Constraining ", nfrag," fragments."
604 write(iout,*) "Length of each fragment, weight and q0:"
606 write (iout,*) "Set of restraints #",iset
608 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
609 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
611 write(iout,*) "constraints between ", npair, "fragments."
612 write(iout,*) "constraint pairs, weights and q0:"
614 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
615 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
617 write(iout,*) "angle constraints within ", nfrag_back,
618 & "backbone fragments."
619 write(iout,*) "fragment, weights:"
621 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
622 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
623 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
626 iset=mod(kolor,nset)+1
629 if(me.eq.king.or..not.out1file)
630 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
633 c------------------------------------------------------------------------------
636 C Read molecular data.
638 implicit real*8 (a-h,o-z)
644 include 'COMMON.IOUNITS'
647 include 'COMMON.INTERACT'
648 include 'COMMON.LOCAL'
649 include 'COMMON.NAMES'
650 include 'COMMON.CHAIN'
651 include 'COMMON.FFIELD'
652 include 'COMMON.SBRIDGE'
653 include 'COMMON.HEADER'
654 include 'COMMON.CONTROL'
655 include 'COMMON.DBASE'
656 include 'COMMON.THREAD'
657 include 'COMMON.CONTACTS'
658 include 'COMMON.TORCNSTR'
659 include 'COMMON.TIME1'
660 include 'COMMON.BOUNDS'
662 include 'COMMON.SETUP'
663 include 'COMMON.SHIELD'
664 character*4 sequence(maxres)
666 double precision x(maxvar)
667 character*256 pdbfile
668 character*400 weightcard
669 character*80 weightcard_t,ucase
670 dimension itype_pdb(maxres)
671 common /pizda/ itype_pdb
672 logical seq_comp,fail
673 double precision energia(0:n_ene)
679 C Read weights of the subsequent energy terms.
680 call card_concat(weightcard)
681 call reada(weightcard,'WLONG',wlong,1.0D0)
682 call reada(weightcard,'WSC',wsc,wlong)
683 call reada(weightcard,'WSCP',wscp,wlong)
684 call reada(weightcard,'WELEC',welec,1.0D0)
685 call reada(weightcard,'WVDWPP',wvdwpp,welec)
686 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
687 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
688 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
689 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
690 call reada(weightcard,'WTURN3',wturn3,1.0D0)
691 call reada(weightcard,'WTURN4',wturn4,1.0D0)
692 call reada(weightcard,'WTURN6',wturn6,1.0D0)
693 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
694 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
695 call reada(weightcard,'WBOND',wbond,1.0D0)
696 call reada(weightcard,'WTOR',wtor,1.0D0)
697 call reada(weightcard,'WTORD',wtor_d,1.0D0)
698 call reada(weightcard,'WANG',wang,1.0D0)
699 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
700 call reada(weightcard,'SCAL14',scal14,0.4D0)
701 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
702 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
703 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
704 call reada(weightcard,'TEMP0',temp0,300.0d0)
705 call reada(weightcard,'WSHIELD',wshield,1.0d0)
706 call reada(weightcard,'WLT',wliptran,0.0D0)
707 call reada(weightcard,'WTUBE',wtube,1.0D0)
708 if (index(weightcard,'SOFT').gt.0) ipot=6
709 C 12/1/95 Added weight for the multi-body term WCORR
710 call reada(weightcard,'WCORRH',wcorr,1.0D0)
711 if (wcorr4.gt.0.0d0) wcorr=wcorr4
731 if(me.eq.king.or..not.out1file)
732 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
733 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
735 10 format (/'Energy-term weights (unscaled):'//
736 & 'WSCC= ',f10.6,' (SC-SC)'/
737 & 'WSCP= ',f10.6,' (SC-p)'/
738 & 'WELEC= ',f10.6,' (p-p electr)'/
739 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
740 & 'WBOND= ',f10.6,' (stretching)'/
741 & 'WANG= ',f10.6,' (bending)'/
742 & 'WSCLOC= ',f10.6,' (SC local)'/
743 & 'WTOR= ',f10.6,' (torsional)'/
744 & 'WTORD= ',f10.6,' (double torsional)'/
745 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
746 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
747 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
748 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
749 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
750 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
751 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
752 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
753 & 'WTURN6= ',f10.6,' (turns, 6th order)')
754 if(me.eq.king.or..not.out1file)then
755 if (wcorr4.gt.0.0d0) then
756 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
757 & 'between contact pairs of peptide groups'
758 write (iout,'(2(a,f5.3/))')
759 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
760 & 'Range of quenching the correlation terms:',2*delt_corr
761 else if (wcorr.gt.0.0d0) then
762 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
763 & 'between contact pairs of peptide groups'
765 write (iout,'(a,f8.3)')
766 & 'Scaling factor of 1,4 SC-p interactions:',scal14
767 write (iout,'(a,f8.3)')
768 & 'General scaling factor of SC-p interactions:',scalscp
770 r0_corr=cutoff_corr-delt_corr
772 aad(i,1)=scalscp*aad(i,1)
773 aad(i,2)=scalscp*aad(i,2)
774 bad(i,1)=scalscp*bad(i,1)
775 bad(i,2)=scalscp*bad(i,2)
777 call rescale_weights(t_bath)
778 if(me.eq.king.or..not.out1file)
779 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
780 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
782 22 format (/'Energy-term weights (scaled):'//
783 & 'WSCC= ',f10.6,' (SC-SC)'/
784 & 'WSCP= ',f10.6,' (SC-p)'/
785 & 'WELEC= ',f10.6,' (p-p electr)'/
786 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
787 & 'WBOND= ',f10.6,' (stretching)'/
788 & 'WANG= ',f10.6,' (bending)'/
789 & 'WSCLOC= ',f10.6,' (SC local)'/
790 & 'WTOR= ',f10.6,' (torsional)'/
791 & 'WTORD= ',f10.6,' (double torsional)'/
792 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
793 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
794 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
795 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
796 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
797 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
798 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
799 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
800 & 'WTURN6= ',f10.6,' (turns, 6th order)')
801 if(me.eq.king.or..not.out1file)
802 & write (iout,*) "Reference temperature for weights calculation:",
804 call reada(weightcard,"D0CM",d0cm,3.78d0)
805 call reada(weightcard,"AKCM",akcm,15.1d0)
806 call reada(weightcard,"AKTH",akth,11.0d0)
807 call reada(weightcard,"AKCT",akct,12.0d0)
808 call reada(weightcard,"V1SS",v1ss,-1.08d0)
809 call reada(weightcard,"V2SS",v2ss,7.61d0)
810 call reada(weightcard,"V3SS",v3ss,13.7d0)
811 call reada(weightcard,"EBR",ebr,-5.50D0)
812 call reada(weightcard,"ATRISS",atriss,0.301D0)
813 call reada(weightcard,"BTRISS",btriss,0.021D0)
814 call reada(weightcard,"CTRISS",ctriss,1.001D0)
815 call reada(weightcard,"DTRISS",dtriss,1.001D0)
816 call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
817 write (iout,*) "ATRISS=", atriss
818 write (iout,*) "BTRISS=", btriss
819 write (iout,*) "CTRISS=", ctriss
820 write (iout,*) "DTRISS=", dtriss
821 if (shield_mode.gt.0) then
823 write (iout,*) "liscale not used in shield mode"
825 write (iout,*) "lipid scaling factor", lipscale
826 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
828 dyn_ss_mask(i)=.false.
832 dyn_ssbond_ij(i,j)=1.0d300
835 call reada(weightcard,"HT",Ht,0.0D0)
837 ss_depth=ebr/wsc-0.25*eps(1,1)
838 Ht=Ht/wsc-0.25*eps(1,1)
846 if (wstrain.ne.0.0) then
847 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
852 write (iout,*) "wshield,", wshield
853 if(me.eq.king.or..not.out1file) then
854 write (iout,*) "Parameters of the SS-bond potential:"
855 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
857 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
858 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
859 write (iout,*)" HT",Ht
860 print *,'indpdb=',indpdb,' pdbref=',pdbref
862 if (indpdb.gt.0 .or. pdbref) then
863 read(inp,'(a)') pdbfile
864 if(me.eq.king.or..not.out1file)
865 & write (iout,'(2a)') 'PDB data will be read from file ',
866 & pdbfile(:ilen(pdbfile))
867 open(ipdbin,file=pdbfile,status='old',err=33)
869 33 write (iout,'(a)') 'Error opening PDB file.'
872 c write (iout,*) 'Begin reading pdb data'
875 c write (iout,*) 'Finished reading pdb data'
877 if(me.eq.king.or..not.out1file)
878 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
879 & ' nstart_sup=',nstart_sup
881 itype_pdb(i)=itype(i)
885 nct=nstart_sup+nsup-1
886 call contact(.false.,ncont_ref,icont_ref,co)
889 if(me.eq.king.or..not.out1file)
890 & write(iout,*)'Adding sidechains'
894 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
897 do while (fail.and.nsi.le.maxsi)
898 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
901 if(fail) write(iout,*)'Adding sidechain failed for res ',
902 & i,' after ',nsi,' trials'
907 if (indpdb.eq.0) then
908 C Read sequence if not taken from the pdb file.
910 c print *,'nres=',nres
911 if (iscode.gt.0) then
912 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
914 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
916 C Convert sequence to numeric code
918 itype(i)=rescode(i,sequence(i),iscode)
920 C Assign initial virtual bond lengths
926 vbld(i+nres)=dsc(iabs(itype(i)))
927 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
928 c write (iout,*) "i",i," itype",itype(i),
929 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
933 c print '(20i4)',(itype(i),i=1,nres)
936 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
938 if (itype(i).eq.ntyp1) then
942 else if (iabs(itype(i+1)).ne.20) then
944 else if (iabs(itype(i)).ne.20) then
951 if(me.eq.king.or..not.out1file)then
952 write (iout,*) "ITEL"
954 write (iout,*) i,itype(i),itel(i)
956 print *,'Call Read_Bridge.'
959 C 8/13/98 Set limits to generating the dihedral angles
964 read (inp,*) ndih_constr
965 if (ndih_constr.gt.0) then
967 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
969 if(me.eq.king.or..not.out1file)then
971 & 'There are',ndih_constr,' constraints on phi angles.'
973 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
978 phi0(i)=deg2rad*phi0(i)
979 drange(i)=deg2rad*drange(i)
981 C if(me.eq.king.or..not.out1file)
982 C & write (iout,*) 'FTORS',ftors
985 phibound(1,ii) = phi0(i)-drange(i)
986 phibound(2,ii) = phi0(i)+drange(i)
989 C first setting the theta boundaries to 0 to pi
990 C this mean that there is no energy penalty for any angle occuring this can be applied
991 C for generate random conformation but is not implemented in this way
996 C begin reading theta constrains this is quartic constrains allowing to
997 C have smooth second derivative
998 if (with_theta_constr) then
999 C with_theta_constr is keyword allowing for occurance of theta constrains
1000 read (inp,*) ntheta_constr
1001 C ntheta_constr is the number of theta constrains
1002 if (ntheta_constr.gt.0) then
1003 C read (inp,*) ftors
1004 read (inp,*) (itheta_constr(i),theta_constr0(i),
1005 & theta_drange(i),for_thet_constr(i),
1006 & i=1,ntheta_constr)
1007 C the above code reads from 1 to ntheta_constr
1008 C itheta_constr(i) residue i for which is theta_constr
1009 C theta_constr0 the global minimum value
1010 C theta_drange is range for which there is no energy penalty
1011 C for_thet_constr is the force constant for quartic energy penalty
1013 if(me.eq.king.or..not.out1file)then
1015 & 'There are',ntheta_constr,' constraints on phi angles.'
1016 do i=1,ntheta_constr
1017 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1019 & for_thet_constr(i)
1022 do i=1,ntheta_constr
1023 theta_constr0(i)=deg2rad*theta_constr0(i)
1024 theta_drange(i)=deg2rad*theta_drange(i)
1026 C if(me.eq.king.or..not.out1file)
1027 C & write (iout,*) 'FTORS',ftors
1028 C do i=1,ntheta_constr
1029 C ii = itheta_constr(i)
1030 C thetabound(1,ii) = phi0(i)-drange(i)
1031 C thetabound(2,ii) = phi0(i)+drange(i)
1033 endif ! ntheta_constr.gt.0
1034 endif! with_theta_constr
1036 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1037 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1040 if (me.eq.king) then
1042 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1044 write (iout,'(a3,i5,2f10.1)')
1045 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1051 cd print *,'NNT=',NNT,' NCT=',NCT
1052 if (itype(1).eq.ntyp1) nnt=2
1053 if (itype(nres).eq.ntyp1) nct=nct-1
1055 if(me.eq.king.or..not.out1file)
1056 & write (iout,'(a,i3)') 'nsup=',nsup
1058 if (nsup.le.(nct-nnt+1)) then
1059 do i=0,nct-nnt+1-nsup
1060 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1066 & 'Error - sequences to be superposed do not match.'
1069 do i=0,nsup-(nct-nnt+1)
1070 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1072 nstart_sup=nstart_sup+i
1078 & 'Error - sequences to be superposed do not match.'
1081 if (nsup.eq.0) nsup=nct-nnt
1082 if (nstart_sup.eq.0) nstart_sup=nnt
1083 if (nstart_seq.eq.0) nstart_seq=nnt
1084 if(me.eq.king.or..not.out1file)
1085 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1086 & ' nstart_seq=',nstart_seq
1088 c--- Zscore rms -------
1089 if (nz_start.eq.0) nz_start=nnt
1090 if (nz_end.eq.0 .and. nsup.gt.0) then
1092 else if (nz_end.eq.0) then
1095 if(me.eq.king.or..not.out1file)then
1096 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1097 write (iout,*) 'IZ_SC=',iz_sc
1099 c----------------------
1102 if (.not.pdbref) then
1103 call read_angles(inp,*38)
1105 38 write (iout,'(a)') 'Error reading reference structure.'
1107 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1108 stop 'Error reading reference structure'
1110 39 call chainbuild_extconf
1112 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1119 cref(j,i,kkk)=c(j,i)
1122 call contact(.true.,ncont_ref,icont_ref,co)
1126 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1128 if (constr_dist.gt.0) call read_dist_constr
1129 write (iout,*) "After read_dist_constr nhpb",nhpb
1130 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1132 if(me.eq.king.or..not.out1file)
1133 & write (iout,*) 'Contact order:',co
1135 if(me.eq.king.or..not.out1file)
1136 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1139 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1141 if(me.eq.king.or..not.out1file)
1142 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1143 & icont_ref(1,i),' ',
1144 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1148 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1149 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1150 & modecalc.ne.10) then
1151 C If input structure hasn't been supplied from the PDB file read or generate
1153 if (iranconf.eq.0 .and. .not. extconf) then
1154 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1155 & write (iout,'(a)') 'Initial geometry will be read in.'
1157 read(inp,'(8f10.5)',end=36,err=36)
1158 & ((c(l,k),l=1,3),k=1,nres),
1159 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1160 write (iout,*) "Exit READ_CART"
1161 write (iout,'(8f10.5)')
1162 & ((c(l,k),l=1,3),k=1,nres),
1163 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1164 call int_from_cart1(.true.)
1165 write (iout,*) "Finish INT_TO_CART"
1168 dc(j,i)=c(j,i+1)-c(j,i)
1169 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1173 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1175 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1176 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1182 call read_angles(inp,*36)
1183 call chainbuild_extconf
1186 36 write (iout,'(a)') 'Error reading angle file.'
1188 call mpi_finalize( MPI_COMM_WORLD,IERR )
1190 stop 'Error reading angle file.'
1192 else if (extconf) then
1193 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1194 & write (iout,'(a)') 'Extended chain initial geometry.'
1196 theta(i)=90d0*deg2rad
1199 phi(i)=180d0*deg2rad
1202 alph(i)=110d0*deg2rad
1205 omeg(i)=-120d0*deg2rad
1206 if (itype(i).le.0) omeg(i)=-omeg(i)
1208 call chainbuild_extconf
1210 if(me.eq.king.or..not.out1file)
1211 & write (iout,'(a)') 'Random-generated initial geometry.'
1215 if (me.eq.king .or. fg_rank.eq.0 .and. (
1216 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1220 call gen_rand_conf(itmp,*30)
1222 30 write (iout,*) 'Failed to generate random conformation',
1223 & ', itrial=',itrial
1224 write (*,*) 'Processor:',me,
1225 & ' Failed to generate random conformation',
1235 write (iout,'(a,i3,a)') 'Processor:',me,
1236 & ' error in generating random conformation.'
1237 write (*,'(a,i3,a)') 'Processor:',me,
1238 & ' error in generating random conformation.'
1241 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1246 & ' error in generating random conformation.'
1251 elseif (modecalc.eq.4) then
1252 read (inp,'(a)') intinname
1253 open (intin,file=intinname,status='old',err=333)
1254 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1255 & write (iout,'(a)') 'intinname',intinname
1256 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1258 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1260 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1262 stop 'Error opening angle file.'
1266 C Generate distance constraints, if the PDB structure is to be regularized.
1267 if (nthread.gt.0) then
1268 call read_threadbase
1271 if (me.eq.king .or. .not. out1file)
1273 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1274 write (iout,'(/a,i3,a)')
1275 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1276 write (iout,'(20i4)') (iss(i),i=1,ns)
1278 write(iout,*)"Running with dynamic disulfide-bond formation"
1280 write (iout,'(/a/)') 'Pre-formed links are:'
1286 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1287 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1293 if (ns.gt.0.and.dyn_ss) then
1297 forcon(i-nss)=forcon(i)
1304 dyn_ss_mask(iss(i))=.true.
1307 if (i2ndstr.gt.0) call secstrp2dihc
1308 c call geom_to_var(nvar,x)
1309 c call etotal(energia(0))
1310 c call enerprint(energia(0))
1311 c call briefout(0,etot)
1313 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1314 cd write (iout,'(a)') 'Variable list:'
1315 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1317 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1318 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1319 & 'Processor',myrank,': end reading molecular data.'
1324 c--------------------------------------------------------------------------
1325 logical function seq_comp(itypea,itypeb,length)
1327 integer length,itypea(length),itypeb(length)
1330 if (itypea(i).ne.itypeb(i)) then
1338 c-----------------------------------------------------------------------------
1339 subroutine read_bridge
1340 C Read information about disulfide bridges.
1341 implicit real*8 (a-h,o-z)
1342 include 'DIMENSIONS'
1346 include 'COMMON.IOUNITS'
1347 include 'COMMON.GEO'
1348 include 'COMMON.VAR'
1349 include 'COMMON.INTERACT'
1350 include 'COMMON.LOCAL'
1351 include 'COMMON.NAMES'
1352 include 'COMMON.CHAIN'
1353 include 'COMMON.FFIELD'
1354 include 'COMMON.SBRIDGE'
1355 include 'COMMON.HEADER'
1356 include 'COMMON.CONTROL'
1357 include 'COMMON.DBASE'
1358 include 'COMMON.THREAD'
1359 include 'COMMON.TIME1'
1360 include 'COMMON.SETUP'
1361 C Read bridging residues.
1362 read (inp,*) ns,(iss(i),i=1,ns)
1364 if(me.eq.king.or..not.out1file)
1365 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1366 C Check whether the specified bridging residues are cystines.
1368 if (itype(iss(i)).ne.1) then
1369 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1370 & 'Do you REALLY think that the residue ',
1371 & restyp(itype(iss(i))),i,
1372 & ' can form a disulfide bridge?!!!'
1373 write (*,'(2a,i3,a)')
1374 & 'Do you REALLY think that the residue ',
1375 & restyp(itype(iss(i))),i,
1376 & ' can form a disulfide bridge?!!!'
1378 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1383 C Read preformed bridges.
1385 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1387 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1390 C Check if the residues involved in bridges are in the specified list of
1391 C bridging residues.
1394 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1395 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1396 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1397 & ' contains residues present in other pairs.'
1398 write (*,'(a,i3,a)') 'Disulfide pair',i,
1399 & ' contains residues present in other pairs.'
1401 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1407 if (ihpb(i).eq.iss(j)) goto 10
1409 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1412 if (jhpb(i).eq.iss(j)) goto 20
1414 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1420 ihpb(i)=ihpb(i)+nres
1421 jhpb(i)=jhpb(i)+nres
1427 c----------------------------------------------------------------------------
1428 subroutine read_x(kanal,*)
1429 implicit real*8 (a-h,o-z)
1430 include 'DIMENSIONS'
1431 include 'COMMON.GEO'
1432 include 'COMMON.VAR'
1433 include 'COMMON.CHAIN'
1434 include 'COMMON.IOUNITS'
1435 include 'COMMON.CONTROL'
1436 include 'COMMON.LOCAL'
1437 include 'COMMON.INTERACT'
1438 c Read coordinates from input
1440 read(kanal,'(8f10.5)',end=10,err=10)
1441 & ((c(l,k),l=1,3),k=1,nres),
1442 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1445 c(j,2*nres)=c(j,nres)
1447 call int_from_cart1(.false.)
1450 dc(j,i)=c(j,i+1)-c(j,i)
1451 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1455 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1457 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1458 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1466 c----------------------------------------------------------------------------
1467 subroutine read_threadbase
1468 implicit real*8 (a-h,o-z)
1469 include 'DIMENSIONS'
1470 include 'COMMON.IOUNITS'
1471 include 'COMMON.GEO'
1472 include 'COMMON.VAR'
1473 include 'COMMON.INTERACT'
1474 include 'COMMON.LOCAL'
1475 include 'COMMON.NAMES'
1476 include 'COMMON.CHAIN'
1477 include 'COMMON.FFIELD'
1478 include 'COMMON.SBRIDGE'
1479 include 'COMMON.HEADER'
1480 include 'COMMON.CONTROL'
1481 include 'COMMON.DBASE'
1482 include 'COMMON.THREAD'
1483 include 'COMMON.TIME1'
1484 C Read pattern database for threading.
1485 read (icbase,*) nseq
1487 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1488 & nres_base(2,i),nres_base(3,i)
1489 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1491 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1492 c & nres_base(2,i),nres_base(3,i)
1493 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1497 if (weidis.eq.0.0D0) weidis=0.1D0
1506 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1507 write (iout,'(a,i5)') 'nexcl: ',nexcl
1508 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1511 c------------------------------------------------------------------------------
1512 subroutine setup_var
1513 implicit real*8 (a-h,o-z)
1514 include 'DIMENSIONS'
1515 include 'COMMON.IOUNITS'
1516 include 'COMMON.GEO'
1517 include 'COMMON.VAR'
1518 include 'COMMON.INTERACT'
1519 include 'COMMON.LOCAL'
1520 include 'COMMON.NAMES'
1521 include 'COMMON.CHAIN'
1522 include 'COMMON.FFIELD'
1523 include 'COMMON.SBRIDGE'
1524 include 'COMMON.HEADER'
1525 include 'COMMON.CONTROL'
1526 include 'COMMON.DBASE'
1527 include 'COMMON.THREAD'
1528 include 'COMMON.TIME1'
1529 C Set up variable list.
1535 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1537 ialph(i,1)=nvar+nside
1541 if (indphi.gt.0) then
1543 else if (indback.gt.0) then
1548 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1551 c----------------------------------------------------------------------------
1552 subroutine gen_dist_constr
1553 C Generate CA distance constraints.
1554 implicit real*8 (a-h,o-z)
1555 include 'DIMENSIONS'
1556 include 'COMMON.IOUNITS'
1557 include 'COMMON.GEO'
1558 include 'COMMON.VAR'
1559 include 'COMMON.INTERACT'
1560 include 'COMMON.LOCAL'
1561 include 'COMMON.NAMES'
1562 include 'COMMON.CHAIN'
1563 include 'COMMON.FFIELD'
1564 include 'COMMON.SBRIDGE'
1565 include 'COMMON.HEADER'
1566 include 'COMMON.CONTROL'
1567 include 'COMMON.DBASE'
1568 include 'COMMON.THREAD'
1569 include 'COMMON.TIME1'
1570 dimension itype_pdb(maxres)
1571 common /pizda/ itype_pdb
1573 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1574 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1575 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1577 if (constr_dist.eq.11) then
1578 do i=nstart_sup,nstart_sup+nsup-1
1579 do j=i+2,nstart_sup+nsup-1
1581 if (distance.le.15.0) then
1583 ihpb(nhpb)=i+nstart_seq-nstart_sup
1584 jhpb(nhpb)=j+nstart_seq-nstart_sup
1585 forcon(nhpb)=sqrt(0.04*distance)
1586 fordepth(nhpb)=sqrt(40.0/distance)
1587 dhpb(nhpb)=distance-0.1d0
1588 dhpb1(nhpb)=distance+0.1d0
1591 if (.not.out1file .or. me.eq.king)
1592 & write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1593 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1595 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ",
1596 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1602 do i=nstart_sup,nstart_sup+nsup-1
1603 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1604 cd & ' seq_pdb', restyp(itype_pdb(i))
1605 do j=i+2,nstart_sup+nsup-1
1607 ihpb(nhpb)=i+nstart_seq-nstart_sup
1608 jhpb(nhpb)=j+nstart_seq-nstart_sup
1610 dhpb(nhpb)=dist(i,j)
1614 cd write (iout,'(a)') 'Distance constraints:'
1619 cd if (ii.gt.nres) then
1624 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1625 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1626 cd & dhpb(i),forcon(i)
1630 c----------------------------------------------------------------------------
1632 implicit real*8 (a-h,o-z)
1633 include 'DIMENSIONS'
1634 include 'COMMON.MAP'
1635 include 'COMMON.IOUNITS'
1636 character*3 angid(4) /'THE','PHI','ALP','OME'/
1637 character*80 mapcard,ucase
1639 read (inp,'(a)') mapcard
1640 mapcard=ucase(mapcard)
1641 if (index(mapcard,'PHI').gt.0) then
1643 else if (index(mapcard,'THE').gt.0) then
1645 else if (index(mapcard,'ALP').gt.0) then
1647 else if (index(mapcard,'OME').gt.0) then
1650 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1651 stop 'Error - illegal variable spec in MAP card.'
1653 call readi (mapcard,'RES1',res1(imap),0)
1654 call readi (mapcard,'RES2',res2(imap),0)
1655 if (res1(imap).eq.0) then
1656 res1(imap)=res2(imap)
1657 else if (res2(imap).eq.0) then
1658 res2(imap)=res1(imap)
1660 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1662 & 'Error - illegal definition of variable group in MAP.'
1663 stop 'Error - illegal definition of variable group in MAP.'
1665 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1666 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1667 call readi(mapcard,'NSTEP',nstep(imap),0)
1668 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1670 & 'Illegal boundary and/or step size specification in MAP.'
1671 stop 'Illegal boundary and/or step size specification in MAP.'
1676 c----------------------------------------------------------------------------
1678 implicit real*8 (a-h,o-z)
1679 include 'DIMENSIONS'
1680 include 'COMMON.IOUNITS'
1681 include 'COMMON.GEO'
1682 include 'COMMON.CSA'
1683 include 'COMMON.BANK'
1684 include 'COMMON.CONTROL'
1686 character*620 mcmcard
1687 call card_concat(mcmcard)
1689 call readi(mcmcard,'NCONF',nconf,50)
1690 call readi(mcmcard,'NADD',nadd,0)
1691 call readi(mcmcard,'JSTART',jstart,1)
1692 call readi(mcmcard,'JEND',jend,1)
1693 call readi(mcmcard,'NSTMAX',nstmax,500000)
1694 call readi(mcmcard,'N0',n0,1)
1695 call readi(mcmcard,'N1',n1,6)
1696 call readi(mcmcard,'N2',n2,4)
1697 call readi(mcmcard,'N3',n3,0)
1698 call readi(mcmcard,'N4',n4,0)
1699 call readi(mcmcard,'N5',n5,0)
1700 call readi(mcmcard,'N6',n6,10)
1701 call readi(mcmcard,'N7',n7,0)
1702 call readi(mcmcard,'N8',n8,0)
1703 call readi(mcmcard,'N9',n9,0)
1704 call readi(mcmcard,'N14',n14,0)
1705 call readi(mcmcard,'N15',n15,0)
1706 call readi(mcmcard,'N16',n16,0)
1707 call readi(mcmcard,'N17',n17,0)
1708 call readi(mcmcard,'N18',n18,0)
1710 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1712 call readi(mcmcard,'NDIFF',ndiff,2)
1713 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1714 call readi(mcmcard,'IS1',is1,1)
1715 call readi(mcmcard,'IS2',is2,8)
1716 call readi(mcmcard,'NRAN0',nran0,4)
1717 call readi(mcmcard,'NRAN1',nran1,2)
1718 call readi(mcmcard,'IRR',irr,1)
1719 call readi(mcmcard,'NSEED',nseed,20)
1720 call readi(mcmcard,'NTOTAL',ntotal,10000)
1721 call reada(mcmcard,'CUT1',cut1,2.0d0)
1722 call reada(mcmcard,'CUT2',cut2,5.0d0)
1723 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1724 call readi(mcmcard,'ICMAX',icmax,3)
1725 call readi(mcmcard,'IRESTART',irestart,0)
1726 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1729 call reada(mcmcard,'DELE',dele,20.0d0)
1730 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1731 call readi(mcmcard,'IREF',iref,0)
1732 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1733 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1734 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1735 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1736 write (iout,*) "NCONF_IN",nconf_in
1739 c----------------------------------------------------------------------------
1740 cfmc subroutine mcmfread
1741 cfmc implicit real*8 (a-h,o-z)
1742 cfmc include 'DIMENSIONS'
1743 cfmc include 'COMMON.MCMF'
1744 cfmc include 'COMMON.IOUNITS'
1745 cfmc include 'COMMON.GEO'
1746 cfmc character*80 ucase
1747 cfmc character*620 mcmcard
1748 cfmc call card_concat(mcmcard)
1750 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1751 cfmc write(iout,*)'MAXRANT=',maxrant
1752 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1753 cfmc write(iout,*)'MAXFAM=',maxfam
1754 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1755 cfmc write(iout,*)'NNET1=',nnet1
1756 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1757 cfmc write(iout,*)'NNET2=',nnet2
1758 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1759 cfmc write(iout,*)'NNET3=',nnet3
1760 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1761 cfmc write(iout,*)'ILASTT=',ilastt
1762 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1763 cfmc write(iout,*)'MAXSTR=',maxstr
1764 cfmc maxstr_f=maxstr/maxfam
1765 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1766 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1767 cfmc write(iout,*)'NMCMF=',nmcmf
1768 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1769 cfmc write(iout,*)'IFOCUS=',ifocus
1770 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1771 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1772 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1773 cfmc write(iout,*)'INTPRT=',intprt
1774 cfmc call readi(mcmcard,'IPRT',iprt,100)
1775 cfmc write(iout,*)'IPRT=',iprt
1776 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1777 cfmc write(iout,*)'IMAXTR=',imaxtr
1778 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1779 cfmc write(iout,*)'MAXEVEN=',maxeven
1780 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1781 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1782 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1783 cfmc write(iout,*)'INIMIN=',inimin
1784 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1785 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1786 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1787 cfmc write(iout,*)'NTHREAD=',nthread
1788 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1789 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1790 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1791 cfmc write(iout,*)'MAXPERT=',maxpert
1792 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1793 cfmc write(iout,*)'IRMSD=',irmsd
1794 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1795 cfmc write(iout,*)'DENEMIN=',denemin
1796 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1797 cfmc write(iout,*)'RCUT1S=',rcut1s
1798 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1799 cfmc write(iout,*)'RCUT1E=',rcut1e
1800 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1801 cfmc write(iout,*)'RCUT2S=',rcut2s
1802 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1803 cfmc write(iout,*)'RCUT2E=',rcut2e
1804 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1805 cfmc write(iout,*)'DPERT1=',d_pert1
1806 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1807 cfmc write(iout,*)'DPERT1A=',d_pert1a
1808 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1809 cfmc write(iout,*)'DPERT2=',d_pert2
1810 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1811 cfmc write(iout,*)'DPERT2A=',d_pert2a
1812 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1813 cfmc write(iout,*)'DPERT2B=',d_pert2b
1814 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1815 cfmc write(iout,*)'DPERT2C=',d_pert2c
1816 cfmc d_pert1=deg2rad*d_pert1
1817 cfmc d_pert1a=deg2rad*d_pert1a
1818 cfmc d_pert2=deg2rad*d_pert2
1819 cfmc d_pert2a=deg2rad*d_pert2a
1820 cfmc d_pert2b=deg2rad*d_pert2b
1821 cfmc d_pert2c=deg2rad*d_pert2c
1822 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1823 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1824 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1825 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1826 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1827 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1828 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1829 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1830 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1831 cfmc write(iout,*)'RCUTINI=',rcutini
1832 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1833 cfmc write(iout,*)'GRAT=',grat
1834 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1835 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1839 c----------------------------------------------------------------------------
1841 implicit real*8 (a-h,o-z)
1842 include 'DIMENSIONS'
1843 include 'COMMON.MCM'
1844 include 'COMMON.MCE'
1845 include 'COMMON.IOUNITS'
1847 character*320 mcmcard
1848 call card_concat(mcmcard)
1849 call readi(mcmcard,'MAXACC',maxacc,100)
1850 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1851 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1852 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1853 call readi(mcmcard,'MAXREPM',maxrepm,200)
1854 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1855 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1856 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1857 call reada(mcmcard,'E_UP',e_up,5.0D0)
1858 call reada(mcmcard,'DELTE',delte,0.1D0)
1859 call readi(mcmcard,'NSWEEP',nsweep,5)
1860 call readi(mcmcard,'NSTEPH',nsteph,0)
1861 call readi(mcmcard,'NSTEPC',nstepc,0)
1862 call reada(mcmcard,'TMIN',tmin,298.0D0)
1863 call reada(mcmcard,'TMAX',tmax,298.0D0)
1864 call readi(mcmcard,'NWINDOW',nwindow,0)
1865 call readi(mcmcard,'PRINT_MC',print_mc,0)
1866 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1867 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1868 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1869 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1870 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1871 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1872 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1873 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1874 if (nwindow.gt.0) then
1875 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1877 winlen(i)=winend(i)-winstart(i)+1
1880 if (tmax.lt.tmin) tmax=tmin
1881 if (tmax.eq.tmin) then
1885 if (nstepc.gt.0 .and. nsteph.gt.0) then
1886 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1887 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1889 C Probabilities of different move types
1890 sumpro_type(0)=0.0D0
1891 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1892 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1893 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1894 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1895 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1896 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1897 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1899 print *,'i',i,' sumprotype',sumpro_type(i)
1900 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1901 print *,'i',i,' sumprotype',sumpro_type(i)
1905 c----------------------------------------------------------------------------
1906 subroutine read_minim
1907 implicit real*8 (a-h,o-z)
1908 include 'DIMENSIONS'
1909 include 'COMMON.MINIM'
1910 include 'COMMON.IOUNITS'
1912 character*320 minimcard
1913 call card_concat(minimcard)
1914 call readi(minimcard,'MAXMIN',maxmin,2000)
1915 call readi(minimcard,'MAXFUN',maxfun,5000)
1916 call readi(minimcard,'MINMIN',minmin,maxmin)
1917 call readi(minimcard,'MINFUN',minfun,maxmin)
1918 call reada(minimcard,'TOLF',tolf,1.0D-2)
1919 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1920 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1921 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1922 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1923 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1924 & 'Options in energy minimization:'
1925 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1926 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1927 & 'MinMin:',MinMin,' MinFun:',MinFun,
1928 & ' TolF:',TolF,' RTolF:',RTolF
1931 c----------------------------------------------------------------------------
1932 subroutine read_angles(kanal,*)
1933 implicit real*8 (a-h,o-z)
1934 include 'DIMENSIONS'
1935 include 'COMMON.GEO'
1936 include 'COMMON.VAR'
1937 include 'COMMON.CHAIN'
1938 include 'COMMON.IOUNITS'
1939 include 'COMMON.CONTROL'
1940 c Read angles from input
1942 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1943 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1944 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1945 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1948 c 9/7/01 avoid 180 deg valence angle
1949 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1951 theta(i)=deg2rad*theta(i)
1952 phi(i)=deg2rad*phi(i)
1953 alph(i)=deg2rad*alph(i)
1954 omeg(i)=deg2rad*omeg(i)
1959 c----------------------------------------------------------------------------
1960 subroutine reada(rekord,lancuch,wartosc,default)
1962 character*(*) rekord,lancuch
1963 double precision wartosc,default
1966 iread=index(rekord,lancuch)
1967 if (iread.eq.0) then
1971 iread=iread+ilen(lancuch)+1
1972 read (rekord(iread:),*,err=10,end=10) wartosc
1977 c----------------------------------------------------------------------------
1978 subroutine readi(rekord,lancuch,wartosc,default)
1980 character*(*) rekord,lancuch
1981 integer wartosc,default
1984 iread=index(rekord,lancuch)
1985 if (iread.eq.0) then
1989 iread=iread+ilen(lancuch)+1
1990 read (rekord(iread:),*,err=10,end=10) wartosc
1995 c----------------------------------------------------------------------------
1996 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1999 integer tablica(dim),default
2000 character*(*) rekord,lancuch
2007 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2008 if (iread.eq.0) return
2009 iread=iread+ilen(lancuch)+1
2010 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2013 c----------------------------------------------------------------------------
2014 subroutine multreada(rekord,lancuch,tablica,dim,default)
2017 double precision tablica(dim),default
2018 character*(*) rekord,lancuch
2025 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2026 if (iread.eq.0) return
2027 iread=iread+ilen(lancuch)+1
2028 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2031 c----------------------------------------------------------------------------
2032 subroutine openunits
2033 implicit real*8 (a-h,o-z)
2034 include 'DIMENSIONS'
2037 character*16 form,nodename
2040 include 'COMMON.SETUP'
2041 include 'COMMON.IOUNITS'
2043 include 'COMMON.CONTROL'
2044 integer lenpre,lenpot,ilen,lentmp
2046 character*3 out1file_text,ucase
2049 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2050 call getenv_loc("PREFIX",prefix)
2052 call getenv_loc("POT",pot)
2053 call getenv_loc("DIRTMP",tmpdir)
2054 call getenv_loc("CURDIR",curdir)
2055 call getenv_loc("OUT1FILE",out1file_text)
2056 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2057 out1file_text=ucase(out1file_text)
2058 if (out1file_text(1:1).eq."Y") then
2061 out1file=fg_rank.gt.0
2066 if (lentmp.gt.0) then
2067 write (*,'(80(1h!))')
2068 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2069 write (*,'(80(1h!))')
2070 write (*,*)"All output files will be on node /tmp directory."
2072 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2073 if (me.eq.king) then
2074 write (*,*) "The master node is ",nodename
2075 else if (fg_rank.eq.0) then
2076 write (*,*) "I am the CG slave node ",nodename
2078 write (*,*) "I am the FG slave node ",nodename
2081 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2082 lenpre = lentmp+lenpre+1
2084 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2085 C Get the names and open the input files
2086 #if defined(WINIFL) || defined(WINPGI)
2087 open(1,file=pref_orig(:ilen(pref_orig))//
2088 & '.inp',status='old',readonly,shared)
2089 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2090 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2091 C Get parameter filenames and open the parameter files.
2092 call getenv_loc('BONDPAR',bondname)
2093 open (ibond,file=bondname,status='old',readonly,shared)
2094 call getenv_loc('THETPAR',thetname)
2095 open (ithep,file=thetname,status='old',readonly,shared)
2096 call getenv_loc('ROTPAR',rotname)
2097 open (irotam,file=rotname,status='old',readonly,shared)
2098 call getenv_loc('TORPAR',torname)
2099 open (itorp,file=torname,status='old',readonly,shared)
2100 call getenv_loc('TORDPAR',tordname)
2101 open (itordp,file=tordname,status='old',readonly,shared)
2102 call getenv_loc('TORKCC',torkccname)
2103 open (itorkcc,file=torkccname,status='old',readonly,shared)
2104 call getenv_loc('THETKCC',thetkccname)
2105 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2106 call getenv_loc('FOURIER',fouriername)
2107 open (ifourier,file=fouriername,status='old',readonly,shared)
2108 call getenv_loc('ELEPAR',elename)
2109 open (ielep,file=elename,status='old',readonly,shared)
2110 call getenv_loc('SIDEPAR',sidename)
2111 open (isidep,file=sidename,status='old',readonly,shared)
2112 call getenv_loc('LIPTRANPAR',liptranname)
2113 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2114 #elif (defined CRAY) || (defined AIX)
2115 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2117 c print *,"Processor",myrank," opened file 1"
2118 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2119 c print *,"Processor",myrank," opened file 9"
2120 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2121 C Get parameter filenames and open the parameter files.
2122 call getenv_loc('BONDPAR',bondname)
2123 open (ibond,file=bondname,status='old',action='read')
2124 c print *,"Processor",myrank," opened file IBOND"
2125 call getenv_loc('THETPAR',thetname)
2126 open (ithep,file=thetname,status='old',action='read')
2127 c print *,"Processor",myrank," opened file ITHEP"
2128 call getenv_loc('ROTPAR',rotname)
2129 open (irotam,file=rotname,status='old',action='read')
2130 c print *,"Processor",myrank," opened file IROTAM"
2131 call getenv_loc('TORPAR',torname)
2132 open (itorp,file=torname,status='old',action='read')
2133 c print *,"Processor",myrank," opened file ITORP"
2134 call getenv_loc('TORDPAR',tordname)
2135 open (itordp,file=tordname,status='old',action='read')
2136 call getenv_loc('TORKCC',torkccname)
2137 open (itorkcc,file=torkccname,status='old',action='read')
2138 call getenv_loc('THETKCC',thetkccname)
2139 open (ithetkcc,file=thetkccname,status='old',action='read')
2140 c print *,"Processor",myrank," opened file ITORDP"
2141 call getenv_loc('SCCORPAR',sccorname)
2142 open (isccor,file=sccorname,status='old',action='read')
2143 c print *,"Processor",myrank," opened file ISCCOR"
2144 call getenv_loc('FOURIER',fouriername)
2145 open (ifourier,file=fouriername,status='old',action='read')
2146 c print *,"Processor",myrank," opened file IFOURIER"
2147 call getenv_loc('ELEPAR',elename)
2148 open (ielep,file=elename,status='old',action='read')
2149 c print *,"Processor",myrank," opened file IELEP"
2150 call getenv_loc('SIDEPAR',sidename)
2151 open (isidep,file=sidename,status='old',action='read')
2152 call getenv_loc('LIPTRANPAR',liptranname)
2153 open (iliptranpar,file=liptranname,status='old',action='read')
2154 c print *,"Processor",myrank," opened file ISIDEP"
2155 c print *,"Processor",myrank," opened parameter files"
2157 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2158 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2159 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2160 C Get parameter filenames and open the parameter files.
2161 call getenv_loc('BONDPAR',bondname)
2162 open (ibond,file=bondname,status='old')
2163 call getenv_loc('THETPAR',thetname)
2164 open (ithep,file=thetname,status='old')
2165 call getenv_loc('ROTPAR',rotname)
2166 open (irotam,file=rotname,status='old')
2167 call getenv_loc('TORPAR',torname)
2168 open (itorp,file=torname,status='old')
2169 call getenv_loc('TORDPAR',tordname)
2170 open (itordp,file=tordname,status='old')
2171 call getenv_loc('TORKCC',torkccname)
2172 open (itorkcc,file=torkccname,status='old')
2173 call getenv_loc('THETKCC',thetkccname)
2174 open (ithetkcc,file=thetkccname,status='old')
2175 call getenv_loc('SCCORPAR',sccorname)
2176 open (isccor,file=sccorname,status='old')
2177 call getenv_loc('FOURIER',fouriername)
2178 open (ifourier,file=fouriername,status='old')
2179 call getenv_loc('ELEPAR',elename)
2180 open (ielep,file=elename,status='old')
2181 call getenv_loc('SIDEPAR',sidename)
2182 open (isidep,file=sidename,status='old')
2183 call getenv_loc('LIPTRANPAR',liptranname)
2184 open (iliptranpar,file=liptranname,status='old')
2186 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2188 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2189 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2190 C Get parameter filenames and open the parameter files.
2191 call getenv_loc('BONDPAR',bondname)
2192 open (ibond,file=bondname,status='old',readonly)
2193 call getenv_loc('THETPAR',thetname)
2194 open (ithep,file=thetname,status='old',readonly)
2195 call getenv_loc('ROTPAR',rotname)
2196 open (irotam,file=rotname,status='old',readonly)
2197 call getenv_loc('TORPAR',torname)
2198 open (itorp,file=torname,status='old',readonly)
2199 call getenv_loc('TORDPAR',tordname)
2200 open (itordp,file=tordname,status='old',readonly)
2201 call getenv_loc('TORKCC',torkccname)
2202 open (itorkcc,file=torkccname,status='old',readonly)
2203 call getenv_loc('THETKCC',thetkccname)
2204 open (ithetkcc,file=thetkccname,status='old',readonly)
2205 call getenv_loc('SCCORPAR',sccorname)
2206 open (isccor,file=sccorname,status='old',readonly)
2208 call getenv_loc('THETPARPDB',thetname_pdb)
2209 print *,"thetname_pdb ",thetname_pdb
2210 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2211 print *,ithep_pdb," opened"
2213 call getenv_loc('FOURIER',fouriername)
2214 open (ifourier,file=fouriername,status='old',readonly)
2215 call getenv_loc('ELEPAR',elename)
2216 open (ielep,file=elename,status='old',readonly)
2217 call getenv_loc('SIDEPAR',sidename)
2218 open (isidep,file=sidename,status='old',readonly)
2219 call getenv_loc('LIPTRANPAR',liptranname)
2220 open (iliptranpar,file=liptranname,status='old',action='read')
2222 call getenv_loc('ROTPARPDB',rotname_pdb)
2223 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2226 call getenv_loc('TUBEPAR',tubename)
2227 open (itube,file=tubename,status='old',readonly)
2230 C 8/9/01 In the newest version SCp interaction constants are read from a file
2231 C Use -DOLDSCP to use hard-coded constants instead.
2233 call getenv_loc('SCPPAR',scpname)
2234 #if defined(WINIFL) || defined(WINPGI)
2235 open (iscpp,file=scpname,status='old',readonly,shared)
2236 #elif (defined CRAY) || (defined AIX)
2237 open (iscpp,file=scpname,status='old',action='read')
2239 open (iscpp,file=scpname,status='old')
2241 open (iscpp,file=scpname,status='old',readonly)
2244 call getenv_loc('PATTERN',patname)
2245 #if defined(WINIFL) || defined(WINPGI)
2246 open (icbase,file=patname,status='old',readonly,shared)
2247 #elif (defined CRAY) || (defined AIX)
2248 open (icbase,file=patname,status='old',action='read')
2250 open (icbase,file=patname,status='old')
2252 open (icbase,file=patname,status='old',readonly)
2255 C Open output file only for CG processes
2256 c print *,"Processor",myrank," fg_rank",fg_rank
2257 if (fg_rank.eq.0) then
2259 if (nodes.eq.1) then
2262 npos = dlog10(dfloat(nodes-1))+1
2264 if (npos.lt.3) npos=3
2265 write (liczba,'(i1)') npos
2266 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2268 write (liczba,form) me
2269 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2270 & liczba(:ilen(liczba))
2271 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2273 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2275 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2276 & liczba(:ilen(liczba))//'.mol2'
2277 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2278 & liczba(:ilen(liczba))//'.stat'
2280 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2281 & //liczba(:ilen(liczba))//'.stat')
2282 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2285 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2286 & liczba(:ilen(liczba))//'.const'
2291 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2292 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2293 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2294 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2295 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2297 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2299 rest2name=prefix(:ilen(prefix))//'.rst'
2301 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2304 #if defined(AIX) || defined(PGI)
2305 if (me.eq.king .or. .not. out1file)
2306 & open(iout,file=outname,status='unknown')
2308 if (fg_rank.gt.0) then
2309 write (liczba,'(i3.3)') myrank/nfgtasks
2310 write (ll,'(bz,i3.3)') fg_rank
2311 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2316 open(igeom,file=intname,status='unknown',position='append')
2317 open(ipdb,file=pdbname,status='unknown')
2318 open(imol2,file=mol2name,status='unknown')
2319 open(istat,file=statname,status='unknown',position='append')
2321 c1out open(iout,file=outname,status='unknown')
2324 if (me.eq.king .or. .not.out1file)
2325 & open(iout,file=outname,status='unknown')
2327 if (fg_rank.gt.0) then
2328 write (liczba,'(i3.3)') myrank/nfgtasks
2329 write (ll,'(bz,i3.3)') fg_rank
2330 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2335 open(igeom,file=intname,status='unknown',access='append')
2336 open(ipdb,file=pdbname,status='unknown')
2337 open(imol2,file=mol2name,status='unknown')
2338 open(istat,file=statname,status='unknown',access='append')
2340 c1out open(iout,file=outname,status='unknown')
2343 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2344 csa_seed=prefix(:lenpre)//'.CSA.seed'
2345 csa_history=prefix(:lenpre)//'.CSA.history'
2346 csa_bank=prefix(:lenpre)//'.CSA.bank'
2347 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2348 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2349 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2350 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2351 csa_int=prefix(:lenpre)//'.int'
2352 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2353 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2354 csa_in=prefix(:lenpre)//'.CSA.in'
2355 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2358 write (iout,'(80(1h-))')
2359 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2360 write (iout,'(80(1h-))')
2361 write (iout,*) "Input file : ",
2362 & pref_orig(:ilen(pref_orig))//'.inp'
2363 write (iout,*) "Output file : ",
2364 & outname(:ilen(outname))
2366 write (iout,*) "Sidechain potential file : ",
2367 & sidename(:ilen(sidename))
2369 write (iout,*) "SCp potential file : ",
2370 & scpname(:ilen(scpname))
2372 write (iout,*) "Electrostatic potential file : ",
2373 & elename(:ilen(elename))
2374 write (iout,*) "Cumulant coefficient file : ",
2375 & fouriername(:ilen(fouriername))
2376 write (iout,*) "Torsional parameter file : ",
2377 & torname(:ilen(torname))
2378 write (iout,*) "Double torsional parameter file : ",
2379 & tordname(:ilen(tordname))
2380 write (iout,*) "SCCOR parameter file : ",
2381 & sccorname(:ilen(sccorname))
2382 write (iout,*) "Bond & inertia constant file : ",
2383 & bondname(:ilen(bondname))
2384 write (iout,*) "Bending parameter file : ",
2385 & thetname(:ilen(thetname))
2386 write (iout,*) "Rotamer parameter file : ",
2387 & rotname(:ilen(rotname))
2388 write (iout,*) "Thetpdb parameter file : ",
2389 & thetname_pdb(:ilen(thetname_pdb))
2390 write (iout,*) "Threading database : ",
2391 & patname(:ilen(patname))
2393 &write (iout,*)" DIRTMP : ",
2395 write (iout,'(80(1h-))')
2399 c----------------------------------------------------------------------------
2400 subroutine card_concat(card)
2401 implicit real*8 (a-h,o-z)
2402 include 'DIMENSIONS'
2403 include 'COMMON.IOUNITS'
2405 character*80 karta,ucase
2407 read (inp,'(a)') karta
2410 do while (karta(80:80).eq.'&')
2411 card=card(:ilen(card)+1)//karta(:79)
2412 read (inp,'(a)') karta
2415 card=card(:ilen(card)+1)//karta
2418 c----------------------------------------------------------------------------------
2420 implicit real*8 (a-h,o-z)
2421 include 'DIMENSIONS'
2422 include 'COMMON.CHAIN'
2423 include 'COMMON.IOUNITS'
2425 open(irest2,file=rest2name,status='unknown')
2426 read(irest2,*) totT,EK,potE,totE,t_bath
2429 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2432 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2435 read (irest2,*) iset
2440 c---------------------------------------------------------------------------------
2441 subroutine read_fragments
2442 implicit real*8 (a-h,o-z)
2443 include 'DIMENSIONS'
2447 include 'COMMON.SETUP'
2448 include 'COMMON.CHAIN'
2449 include 'COMMON.IOUNITS'
2451 include 'COMMON.CONTROL'
2452 read(inp,*) nset,nfrag,npair,nfrag_back
2453 if(me.eq.king.or..not.out1file)
2454 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2455 & " nfrag_back",nfrag_back
2457 read(inp,*) mset(iset)
2459 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2461 if(me.eq.king.or..not.out1file)
2462 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2463 & ifrag(2,i,iset), qinfrag(i,iset)
2466 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2468 if(me.eq.king.or..not.out1file)
2469 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2470 & ipair(2,i,iset), qinpair(i,iset)
2473 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2474 & wfrag_back(3,i,iset),
2475 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2476 if(me.eq.king.or..not.out1file)
2477 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2478 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2483 C---------------------------------------------------------------------------
2484 subroutine read_afminp
2485 implicit real*8 (a-h,o-z)
2486 include 'DIMENSIONS'
2490 include 'COMMON.SETUP'
2491 include 'COMMON.CONTROL'
2492 include 'COMMON.CHAIN'
2493 include 'COMMON.IOUNITS'
2494 include 'COMMON.SBRIDGE'
2495 character*320 afmcard
2497 call card_concat(afmcard)
2498 call readi(afmcard,"BEG",afmbeg,0)
2499 call readi(afmcard,"END",afmend,0)
2500 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2501 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2502 print *,'FORCE=' ,forceAFMconst
2503 CCCC NOW PROPERTIES FOR AFM
2506 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2508 distafminit=dsqrt(distafminit)
2509 print *,'initdist',distafminit
2513 c-------------------------------------------------------------------------------
2514 subroutine read_dist_constr
2515 implicit real*8 (a-h,o-z)
2516 include 'DIMENSIONS'
2520 include 'COMMON.SETUP'
2521 include 'COMMON.CONTROL'
2522 include 'COMMON.CHAIN'
2523 include 'COMMON.IOUNITS'
2524 include 'COMMON.SBRIDGE'
2525 integer ifrag_(2,100),ipair_(2,100)
2526 double precision wfrag_(100),wpair_(100)
2527 character*500 controlcard
2529 write (iout,*) "Calling read_dist_constr"
2530 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2532 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
2533 call gen_dist_constr
2536 call card_concat(controlcard)
2537 call readi(controlcard,"NFRAG",nfrag_,0)
2538 call readi(controlcard,"NPAIR",npair_,0)
2539 call readi(controlcard,"NDIST",ndist_,0)
2540 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2541 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2542 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2543 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2544 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2545 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2546 c write (iout,*) "IFRAG"
2548 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2550 c write (iout,*) "IPAIR"
2552 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2556 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2557 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2558 & ifrag_(2,i)=nstart_sup+nsup-1
2559 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2561 if (wfrag_(i).gt.0.0d0) then
2562 do j=ifrag_(1,i),ifrag_(2,i)-1
2563 do k=j+1,ifrag_(2,i)
2564 c write (iout,*) "j",j," k",k
2566 if (constr_dist.eq.1) then
2571 forcon(nhpb)=wfrag_(i)
2572 else if (constr_dist.eq.2) then
2573 if (ddjk.le.dist_cut) then
2578 forcon(nhpb)=wfrag_(i)
2585 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2588 if (.not.out1file .or. me.eq.king)
2589 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2590 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2592 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2593 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2600 if (wpair_(i).gt.0.0d0) then
2608 do j=ifrag_(1,ii),ifrag_(2,ii)
2609 do k=ifrag_(1,jj),ifrag_(2,jj)
2613 forcon(nhpb)=wpair_(i)
2614 dhpb(nhpb)=dist(j,k)
2616 if (.not.out1file .or. me.eq.king)
2617 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2618 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2620 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2621 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2629 if (constr_dist.eq.11) then
2630 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2631 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2632 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2635 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2636 & ibecarb(i),forcon(nhpb+1)
2638 if (forcon(nhpb+1).gt.0.0d0) then
2640 if (ibecarb(i).gt.0) then
2641 ihpb(i)=ihpb(i)+nres
2642 jhpb(i)=jhpb(i)+nres
2644 if (dhpb(nhpb).eq.0.0d0)
2645 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2647 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2648 C if (forcon(nhpb+1).gt.0.0d0) then
2650 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2652 if (.not.out1file .or. me.eq.king)
2653 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2654 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2656 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2657 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2665 c-------------------------------------------------------------------------------
2667 subroutine flush(iu)
2672 subroutine flush(iu)
2677 c------------------------------------------------------------------------------
2678 subroutine copy_to_tmp(source)
2679 include "DIMENSIONS"
2680 include "COMMON.IOUNITS"
2681 character*(*) source
2682 character* 256 tmpfile
2686 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2687 inquire(file=tmpfile,exist=ex)
2689 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2690 & " to temporary directory..."
2691 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2692 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2696 c------------------------------------------------------------------------------
2697 subroutine move_from_tmp(source)
2698 include "DIMENSIONS"
2699 include "COMMON.IOUNITS"
2700 character*(*) source
2703 write (*,*) "Moving ",source(:ilen(source)),
2704 & " from temporary directory to working directory"
2705 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2706 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2709 c------------------------------------------------------------------------------
2710 subroutine random_init(seed)
2712 C Initialize random number generator
2714 implicit real*8 (a-h,o-z)
2715 include 'DIMENSIONS'
2718 logical OKRandom, prng_restart
2720 integer iseed_array(4)
2722 include 'COMMON.IOUNITS'
2723 include 'COMMON.TIME1'
2724 include 'COMMON.THREAD'
2725 include 'COMMON.SBRIDGE'
2726 include 'COMMON.CONTROL'
2727 include 'COMMON.MCM'
2728 include 'COMMON.MAP'
2729 include 'COMMON.HEADER'
2730 include 'COMMON.CSA'
2731 include 'COMMON.CHAIN'
2732 include 'COMMON.MUCA'
2734 include 'COMMON.FFIELD'
2735 include 'COMMON.SETUP'
2736 iseed=-dint(dabs(seed))
2737 if (iseed.eq.0) then
2738 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2739 & 'Random seed undefined. The program will stop.'
2740 write (*,'(/80(1h*)/20x,a/80(1h*))')
2741 & 'Random seed undefined. The program will stop.'
2743 call mpi_finalize(mpi_comm_world,ierr)
2745 stop 'Bad random seed.'
2748 if (fg_rank.eq.0) then
2752 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2753 OKRandom = prng_restart(me,iseed)
2756 tmp=65536.0d0**(4-i)
2757 iseed_array(i) = dint(seed/tmp)
2758 seed=seed-iseed_array(i)*tmp
2761 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2762 & (iseed_array(i),i=1,4)
2763 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2764 & (iseed_array(i),i=1,4)
2765 OKRandom = prng_restart(me,iseed_array)
2768 c r1 = prng_next(me)
2769 r1=ran_number(0.0D0,1.0D0)
2771 & write (iout,*) 'ran_num',r1
2772 if (r1.lt.0.0d0) OKRandom=.false.
2774 if (.not.OKRandom) then
2775 write (iout,*) 'PRNG IS NOT WORKING!!!'
2776 print *,'PRNG IS NOT WORKING!!!'
2779 call mpi_abort(mpi_comm_world,error_msg,ierr)
2782 write (iout,*) 'too many processors for parallel prng'
2783 write (*,*) 'too many processors for parallel prng'
2791 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)