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,'IPRINT',iprint,0)
151 C SHIELD keyword sets if the shielding effect of side-chains is used
152 C 0 denots no shielding is used all peptide are equally despite the
153 C solvent accesible area
154 C 1 the newly introduced function
155 C 2 reseved for further possible developement
156 call readi(controlcard,'SHIELD',shield_mode,0)
157 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
158 write(iout,*) "shield_mode",shield_mode
160 call readi(controlcard,'TORMODE',tor_mode,0)
161 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
162 write(iout,*) "torsional and valence angle mode",tor_mode
163 call readi(controlcard,'MAXGEN',maxgen,10000)
164 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
165 call readi(controlcard,"KDIAG",kdiag,0)
166 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
167 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
168 & write (iout,*) "RESCALE_MODE",rescale_mode
169 split_ene=index(controlcard,'SPLIT_ENE').gt.0
170 if (index(controlcard,'REGULAR').gt.0.0D0) then
171 call reada(controlcard,'WEIDIS',weidis,0.1D0)
175 if (index(controlcard,'CHECKGRAD').gt.0) then
177 if (index(controlcard,'CART').gt.0) then
179 elseif (index(controlcard,'CARINT').gt.0) then
184 call reada(controlcard,'DELTA',aincr,1.0d-7)
185 elseif (index(controlcard,'THREAD').gt.0) then
187 call readi(controlcard,'THREAD',nthread,0)
188 if (nthread.gt.0) then
189 call reada(controlcard,'WEIDIS',weidis,0.1D0)
192 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
193 stop 'Error termination in Read_Control.'
195 else if (index(controlcard,'MCMA').gt.0) then
197 else if (index(controlcard,'MCEE').gt.0) then
199 else if (index(controlcard,'MULTCONF').gt.0) then
201 else if (index(controlcard,'MAP').gt.0) then
203 call readi(controlcard,'MAP',nmap,0)
204 else if (index(controlcard,'CSA').gt.0) then
206 crc else if (index(controlcard,'ZSCORE').gt.0) then
208 crc ZSCORE is rm from UNRES, modecalc=9 is available
211 cfcm else if (index(controlcard,'MCMF').gt.0) then
213 else if (index(controlcard,'SOFTREG').gt.0) then
215 else if (index(controlcard,'CHECK_BOND').gt.0) then
217 else if (index(controlcard,'TEST').gt.0) then
219 else if (index(controlcard,'MD').gt.0) then
221 else if (index(controlcard,'RE ').gt.0) then
225 lmuca=index(controlcard,'MUCA').gt.0
226 call readi(controlcard,'MUCADYN',mucadyn,0)
227 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
228 if (lmuca .and. (me.eq.king .or. .not.out1file ))
230 write (iout,*) 'MUCADYN=',mucadyn
231 write (iout,*) 'MUCASMOOTH=',muca_smooth
234 iscode=index(controlcard,'ONE_LETTER')
235 indphi=index(controlcard,'PHI')
236 indback=index(controlcard,'BACK')
237 iranconf=index(controlcard,'RAND_CONF')
238 i2ndstr=index(controlcard,'USE_SEC_PRED')
239 gradout=index(controlcard,'GRADOUT').gt.0
240 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
241 C DISTCHAINMAX become obsolete for periodic boundry condition
242 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
243 C Reading the dimensions of box in x,y,z coordinates
244 call reada(controlcard,'BOXX',boxxsize,100.0d0)
245 call reada(controlcard,'BOXY',boxysize,100.0d0)
246 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
247 c Cutoff range for interactions
248 call reada(controlcard,"R_CUT",r_cut,15.0d0)
249 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
250 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
251 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
252 if (lipthick.gt.0.0d0) then
253 bordliptop=(boxzsize+lipthick)/2.0
254 bordlipbot=bordliptop-lipthick
256 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
257 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
258 buflipbot=bordlipbot+lipbufthick
259 bufliptop=bordliptop-lipbufthick
260 if ((lipbufthick*2.0d0).gt.lipthick)
261 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
263 write(iout,*) "bordliptop=",bordliptop
264 write(iout,*) "bordlipbot=",bordlipbot
265 write(iout,*) "bufliptop=",bufliptop
266 write(iout,*) "buflipbot=",buflipbot
267 write (iout,*) "SHIELD MODE",shield_mode
268 if (TUBElog.gt.0) then
269 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
270 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
271 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
272 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
273 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
274 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
275 buftubebot=bordtubebot+tubebufthick
276 buftubetop=bordtubetop-tubebufthick
278 if (shield_mode.gt.0) then
280 C VSolvSphere the volume of solving sphere
282 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
283 C there will be no distinction between proline peptide group and normal peptide
284 C group in case of shielding parameters
285 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
286 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
287 write (iout,*) VSolvSphere,VSolvSphere_div
288 C long axis of side chain
290 long_r_sidechain(i)=vbldsc0(1,i)
291 short_r_sidechain(i)=sigma0(i)
295 if (me.eq.king .or. .not.out1file )
296 & write (iout,*) "DISTCHAINMAX",distchainmax
298 if(me.eq.king.or..not.out1file)
299 & write (iout,'(2a)') diagmeth(kdiag),
300 & ' routine used to diagonalize matrices.'
303 c--------------------------------------------------------------------------
304 subroutine read_REMDpar
308 implicit real*8 (a-h,o-z)
310 include 'COMMON.IOUNITS'
311 include 'COMMON.TIME1'
314 include 'COMMON.LANGEVIN'
316 include 'COMMON.LANGEVIN.lang0'
318 include 'COMMON.INTERACT'
319 include 'COMMON.NAMES'
321 include 'COMMON.REMD'
322 include 'COMMON.CONTROL'
323 include 'COMMON.SETUP'
325 character*320 controlcard
326 character*3200 controlcard1
327 integer iremd_m_total
329 if(me.eq.king.or..not.out1file)
330 & write (iout,*) "REMD setup"
332 call card_concat(controlcard)
333 call readi(controlcard,"NREP",nrep,3)
334 call readi(controlcard,"NSTEX",nstex,1000)
335 call reada(controlcard,"RETMIN",retmin,10.0d0)
336 call reada(controlcard,"RETMAX",retmax,1000.0d0)
337 mremdsync=(index(controlcard,'SYNC').gt.0)
338 call readi(controlcard,"NSYN",i_sync_step,100)
339 restart1file=(index(controlcard,'REST1FILE').gt.0)
340 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
341 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
342 if(max_cache_traj_use.gt.max_cache_traj)
343 & max_cache_traj_use=max_cache_traj
344 if(me.eq.king.or..not.out1file) then
345 cd if (traj1file) then
346 crc caching is in testing - NTWX is not ignored
347 cd write (iout,*) "NTWX value is ignored"
348 cd write (iout,*) " trajectory is stored to one file by master"
349 cd write (iout,*) " before exchange at NSTEX intervals"
351 write (iout,*) "NREP= ",nrep
352 write (iout,*) "NSTEX= ",nstex
353 write (iout,*) "SYNC= ",mremdsync
354 write (iout,*) "NSYN= ",i_sync_step
355 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
358 if (index(controlcard,'TLIST').gt.0) then
360 call card_concat(controlcard1)
361 read(controlcard1,*) (remd_t(i),i=1,nrep)
362 if(me.eq.king.or..not.out1file)
363 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
366 if (index(controlcard,'MLIST').gt.0) then
368 call card_concat(controlcard1)
369 read(controlcard1,*) (remd_m(i),i=1,nrep)
370 if(me.eq.king.or..not.out1file) then
371 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
374 iremd_m_total=iremd_m_total+remd_m(i)
376 write (iout,*) 'Total number of replicas ',iremd_m_total
379 if(me.eq.king.or..not.out1file)
380 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
383 c--------------------------------------------------------------------------
384 subroutine read_MDpar
388 implicit real*8 (a-h,o-z)
390 include 'COMMON.IOUNITS'
391 include 'COMMON.TIME1'
394 include 'COMMON.LANGEVIN'
396 include 'COMMON.LANGEVIN.lang0'
398 include 'COMMON.INTERACT'
399 include 'COMMON.NAMES'
401 include 'COMMON.SETUP'
402 include 'COMMON.CONTROL'
403 include 'COMMON.SPLITELE'
405 character*320 controlcard
407 call card_concat(controlcard)
408 call readi(controlcard,"NSTEP",n_timestep,1000000)
409 call readi(controlcard,"NTWE",ntwe,100)
410 call readi(controlcard,"NTWX",ntwx,1000)
411 call reada(controlcard,"DT",d_time,1.0d-1)
412 call reada(controlcard,"DVMAX",dvmax,2.0d1)
413 call reada(controlcard,"DAMAX",damax,1.0d1)
414 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
415 call readi(controlcard,"LANG",lang,0)
416 RESPA = index(controlcard,"RESPA") .gt. 0
417 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
418 ntime_split0=ntime_split
419 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
420 ntime_split0=ntime_split
421 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
422 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
423 rest = index(controlcard,"REST").gt.0
424 tbf = index(controlcard,"TBF").gt.0
425 usampl = index(controlcard,"USAMPL").gt.0
427 mdpdb = index(controlcard,"MDPDB").gt.0
428 call reada(controlcard,"T_BATH",t_bath,300.0d0)
429 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
430 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
431 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
432 if (count_reset_moment.eq.0) count_reset_moment=1000000000
433 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
434 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
435 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
436 if (count_reset_vel.eq.0) count_reset_vel=1000000000
437 large = index(controlcard,"LARGE").gt.0
438 print_compon = index(controlcard,"PRINT_COMPON").gt.0
439 rattle = index(controlcard,"RATTLE").gt.0
440 c if performing umbrella sampling, fragments constrained are read from the fragment file
446 if(me.eq.king.or..not.out1file) then
448 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
450 write (iout,'(a)') "The units are:"
451 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
452 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
453 & " acceleration: angstrom/(48.9 fs)**2"
454 write (iout,'(a)') "energy: kcal/mol, temperature: K"
456 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
457 write (iout,'(a60,f10.5,a)')
458 & "Initial time step of numerical integration:",d_time,
460 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
462 write (iout,'(2a,i4,a)')
463 & "A-MTS algorithm used; initial time step for fast-varying",
464 & " short-range forces split into",ntime_split," steps."
465 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
466 & r_cut," lambda",rlamb
468 write (iout,'(2a,f10.5)')
469 & "Maximum acceleration threshold to reduce the time step",
470 & "/increase split number:",damax
471 write (iout,'(2a,f10.5)')
472 & "Maximum predicted energy drift to reduce the timestep",
473 & "/increase split number:",edriftmax
474 write (iout,'(a60,f10.5)')
475 & "Maximum velocity threshold to reduce velocities:",dvmax
476 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
477 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
478 if (rattle) write (iout,'(a60)')
479 & "Rattle algorithm used to constrain the virtual bonds"
483 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
484 call reada(controlcard,"RWAT",rwat,1.4d0)
485 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
486 surfarea=index(controlcard,"SURFAREA").gt.0
487 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
488 if(me.eq.king.or..not.out1file)then
489 write (iout,'(/a,$)') "Langevin dynamics calculation"
492 & " with direct integration of Langevin equations"
493 else if (lang.eq.2) then
494 write (iout,'(a/)') " with TINKER stochasic MD integrator"
495 else if (lang.eq.3) then
496 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
497 else if (lang.eq.4) then
498 write (iout,'(a/)') " in overdamped mode"
500 write (iout,'(//a,i5)')
501 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
504 write (iout,'(a60,f10.5)') "Temperature:",t_bath
505 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
506 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
507 write (iout,'(a60,f10.5)')
508 & "Scaling factor of the friction forces:",scal_fric
509 if (surfarea) write (iout,'(2a,i10,a)')
510 & "Friction coefficients will be scaled by solvent-accessible",
511 & " surface area every",reset_fricmat," steps."
513 c Calculate friction coefficients and bounds of stochastic forces
514 eta=6*pi*cPoise*etawat
515 if(me.eq.king.or..not.out1file)
516 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
518 gamp=scal_fric*(pstok+rwat)*eta
519 stdfp=dsqrt(2*Rb*t_bath/d_time)
521 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
522 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
524 if(me.eq.king.or..not.out1file)then
525 write (iout,'(/2a/)')
526 & "Radii of site types and friction coefficients and std's of",
527 & " stochastic forces of fully exposed sites"
528 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
530 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
531 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
535 if(me.eq.king.or..not.out1file)then
536 write (iout,'(a)') "Berendsen bath calculation"
537 write (iout,'(a60,f10.5)') "Temperature:",t_bath
538 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
540 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
541 & count_reset_moment," steps"
543 & write (iout,'(a,i10,a)')
544 & "Velocities will be reset at random every",count_reset_vel,
548 if(me.eq.king.or..not.out1file)
549 & write (iout,'(a31)') "Microcanonical mode calculation"
551 if(me.eq.king.or..not.out1file)then
552 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
554 write(iout,*) "MD running with constraints."
555 write(iout,*) "Equilibration time ", eq_time, " mtus."
556 write(iout,*) "Constraining ", nfrag," fragments."
557 write(iout,*) "Length of each fragment, weight and q0:"
559 write (iout,*) "Set of restraints #",iset
561 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
562 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
564 write(iout,*) "constraints between ", npair, "fragments."
565 write(iout,*) "constraint pairs, weights and q0:"
567 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
568 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
570 write(iout,*) "angle constraints within ", nfrag_back,
571 & "backbone fragments."
572 write(iout,*) "fragment, weights:"
574 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
575 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
576 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
579 iset=mod(kolor,nset)+1
582 if(me.eq.king.or..not.out1file)
583 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
586 c------------------------------------------------------------------------------
589 C Read molecular data.
591 implicit real*8 (a-h,o-z)
597 include 'COMMON.IOUNITS'
600 include 'COMMON.INTERACT'
601 include 'COMMON.LOCAL'
602 include 'COMMON.NAMES'
603 include 'COMMON.CHAIN'
604 include 'COMMON.FFIELD'
605 include 'COMMON.SBRIDGE'
606 include 'COMMON.HEADER'
607 include 'COMMON.CONTROL'
608 include 'COMMON.DBASE'
609 include 'COMMON.THREAD'
610 include 'COMMON.CONTACTS'
611 include 'COMMON.TORCNSTR'
612 include 'COMMON.TIME1'
613 include 'COMMON.BOUNDS'
615 include 'COMMON.SETUP'
616 include 'COMMON.SHIELD'
617 character*4 sequence(maxres)
619 double precision x(maxvar)
620 character*256 pdbfile
621 character*400 weightcard
622 character*80 weightcard_t,ucase
623 dimension itype_pdb(maxres)
624 common /pizda/ itype_pdb
625 logical seq_comp,fail
626 double precision energia(0:n_ene)
632 C Read weights of the subsequent energy terms.
633 call card_concat(weightcard)
634 call reada(weightcard,'WLONG',wlong,1.0D0)
635 call reada(weightcard,'WSC',wsc,wlong)
636 call reada(weightcard,'WSCP',wscp,wlong)
637 call reada(weightcard,'WELEC',welec,1.0D0)
638 call reada(weightcard,'WVDWPP',wvdwpp,welec)
639 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
640 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
641 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
642 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
643 call reada(weightcard,'WTURN3',wturn3,1.0D0)
644 call reada(weightcard,'WTURN4',wturn4,1.0D0)
645 call reada(weightcard,'WTURN6',wturn6,1.0D0)
646 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
647 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
648 call reada(weightcard,'WBOND',wbond,1.0D0)
649 call reada(weightcard,'WTOR',wtor,1.0D0)
650 call reada(weightcard,'WTORD',wtor_d,1.0D0)
651 call reada(weightcard,'WANG',wang,1.0D0)
652 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
653 call reada(weightcard,'SCAL14',scal14,0.4D0)
654 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
655 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
656 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
657 call reada(weightcard,'TEMP0',temp0,300.0d0)
658 call reada(weightcard,'WSHIELD',wshield,1.0d0)
659 call reada(weightcard,'WLT',wliptran,0.0D0)
660 call reada(weightcard,'WTUBE',wtube,1.0D0)
661 if (index(weightcard,'SOFT').gt.0) ipot=6
662 C 12/1/95 Added weight for the multi-body term WCORR
663 call reada(weightcard,'WCORRH',wcorr,1.0D0)
664 if (wcorr4.gt.0.0d0) wcorr=wcorr4
684 if(me.eq.king.or..not.out1file)
685 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
686 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
688 10 format (/'Energy-term weights (unscaled):'//
689 & 'WSCC= ',f10.6,' (SC-SC)'/
690 & 'WSCP= ',f10.6,' (SC-p)'/
691 & 'WELEC= ',f10.6,' (p-p electr)'/
692 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
693 & 'WBOND= ',f10.6,' (stretching)'/
694 & 'WANG= ',f10.6,' (bending)'/
695 & 'WSCLOC= ',f10.6,' (SC local)'/
696 & 'WTOR= ',f10.6,' (torsional)'/
697 & 'WTORD= ',f10.6,' (double torsional)'/
698 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
699 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
700 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
701 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
702 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
703 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
704 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
705 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
706 & 'WTURN6= ',f10.6,' (turns, 6th order)')
707 if(me.eq.king.or..not.out1file)then
708 if (wcorr4.gt.0.0d0) then
709 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
710 & 'between contact pairs of peptide groups'
711 write (iout,'(2(a,f5.3/))')
712 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
713 & 'Range of quenching the correlation terms:',2*delt_corr
714 else if (wcorr.gt.0.0d0) then
715 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
716 & 'between contact pairs of peptide groups'
718 write (iout,'(a,f8.3)')
719 & 'Scaling factor of 1,4 SC-p interactions:',scal14
720 write (iout,'(a,f8.3)')
721 & 'General scaling factor of SC-p interactions:',scalscp
723 r0_corr=cutoff_corr-delt_corr
725 aad(i,1)=scalscp*aad(i,1)
726 aad(i,2)=scalscp*aad(i,2)
727 bad(i,1)=scalscp*bad(i,1)
728 bad(i,2)=scalscp*bad(i,2)
730 call rescale_weights(t_bath)
731 if(me.eq.king.or..not.out1file)
732 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
733 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
735 22 format (/'Energy-term weights (scaled):'//
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 correlatkion)'/
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)
755 & write (iout,*) "Reference temperature for weights calculation:",
757 call reada(weightcard,"D0CM",d0cm,3.78d0)
758 call reada(weightcard,"AKCM",akcm,15.1d0)
759 call reada(weightcard,"AKTH",akth,11.0d0)
760 call reada(weightcard,"AKCT",akct,12.0d0)
761 call reada(weightcard,"V1SS",v1ss,-1.08d0)
762 call reada(weightcard,"V2SS",v2ss,7.61d0)
763 call reada(weightcard,"V3SS",v3ss,13.7d0)
764 call reada(weightcard,"EBR",ebr,-5.50D0)
765 call reada(weightcard,"ATRISS",atriss,0.301D0)
766 call reada(weightcard,"BTRISS",btriss,0.021D0)
767 call reada(weightcard,"CTRISS",ctriss,1.001D0)
768 call reada(weightcard,"DTRISS",dtriss,1.001D0)
769 write (iout,*) "ATRISS=", atriss
770 write (iout,*) "BTRISS=", btriss
771 write (iout,*) "CTRISS=", ctriss
772 write (iout,*) "DTRISS=", dtriss
773 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
775 dyn_ss_mask(i)=.false.
779 dyn_ssbond_ij(i,j)=1.0d300
782 call reada(weightcard,"HT",Ht,0.0D0)
784 ss_depth=ebr/wsc-0.25*eps(1,1)
785 Ht=Ht/wsc-0.25*eps(1,1)
786 akcm=akcm*wstrain/wsc
787 akth=akth*wstrain/wsc
788 akct=akct*wstrain/wsc
789 v1ss=v1ss*wstrain/wsc
790 v2ss=v2ss*wstrain/wsc
791 v3ss=v3ss*wstrain/wsc
793 if (wstrain.ne.0.0) then
794 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
799 write (iout,*) "wshield,", wshield
800 if(me.eq.king.or..not.out1file) then
801 write (iout,*) "Parameters of the SS-bond potential:"
802 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
804 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
805 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
806 write (iout,*)" HT",Ht
807 print *,'indpdb=',indpdb,' pdbref=',pdbref
809 if (indpdb.gt.0 .or. pdbref) then
810 read(inp,'(a)') pdbfile
811 if(me.eq.king.or..not.out1file)
812 & write (iout,'(2a)') 'PDB data will be read from file ',
813 & pdbfile(:ilen(pdbfile))
814 open(ipdbin,file=pdbfile,status='old',err=33)
816 33 write (iout,'(a)') 'Error opening PDB file.'
819 c write (iout,*) 'Begin reading pdb data'
822 c write (iout,*) 'Finished reading pdb data'
824 if(me.eq.king.or..not.out1file)
825 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
826 & ' nstart_sup=',nstart_sup
828 itype_pdb(i)=itype(i)
832 nct=nstart_sup+nsup-1
833 call contact(.false.,ncont_ref,icont_ref,co)
836 if(me.eq.king.or..not.out1file)
837 & write(iout,*)'Adding sidechains'
841 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
844 do while (fail.and.nsi.le.maxsi)
845 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
848 if(fail) write(iout,*)'Adding sidechain failed for res ',
849 & i,' after ',nsi,' trials'
854 if (indpdb.eq.0) then
855 C Read sequence if not taken from the pdb file.
857 c print *,'nres=',nres
858 if (iscode.gt.0) then
859 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
861 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
863 C Convert sequence to numeric code
865 itype(i)=rescode(i,sequence(i),iscode)
867 C Assign initial virtual bond lengths
873 vbld(i+nres)=dsc(iabs(itype(i)))
874 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
875 c write (iout,*) "i",i," itype",itype(i),
876 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
880 c print '(20i4)',(itype(i),i=1,nres)
883 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
885 if (itype(i).eq.ntyp1) then
889 else if (iabs(itype(i+1)).ne.20) then
891 else if (iabs(itype(i)).ne.20) then
898 if(me.eq.king.or..not.out1file)then
899 write (iout,*) "ITEL"
901 write (iout,*) i,itype(i),itel(i)
903 print *,'Call Read_Bridge.'
906 C 8/13/98 Set limits to generating the dihedral angles
911 read (inp,*) ndih_constr
912 if (ndih_constr.gt.0) then
914 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
916 if(me.eq.king.or..not.out1file)then
918 & 'There are',ndih_constr,' constraints on phi angles.'
920 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
925 phi0(i)=deg2rad*phi0(i)
926 drange(i)=deg2rad*drange(i)
928 C if(me.eq.king.or..not.out1file)
929 C & write (iout,*) 'FTORS',ftors
932 phibound(1,ii) = phi0(i)-drange(i)
933 phibound(2,ii) = phi0(i)+drange(i)
936 C first setting the theta boundaries to 0 to pi
937 C this mean that there is no energy penalty for any angle occuring this can be applied
938 C for generate random conformation but is not implemented in this way
943 C begin reading theta constrains this is quartic constrains allowing to
944 C have smooth second derivative
945 if (with_theta_constr) then
946 C with_theta_constr is keyword allowing for occurance of theta constrains
947 read (inp,*) ntheta_constr
948 C ntheta_constr is the number of theta constrains
949 if (ntheta_constr.gt.0) then
951 read (inp,*) (itheta_constr(i),theta_constr0(i),
952 & theta_drange(i),for_thet_constr(i),
954 C the above code reads from 1 to ntheta_constr
955 C itheta_constr(i) residue i for which is theta_constr
956 C theta_constr0 the global minimum value
957 C theta_drange is range for which there is no energy penalty
958 C for_thet_constr is the force constant for quartic energy penalty
960 if(me.eq.king.or..not.out1file)then
962 & 'There are',ntheta_constr,' constraints on phi angles.'
964 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
970 theta_constr0(i)=deg2rad*theta_constr0(i)
971 theta_drange(i)=deg2rad*theta_drange(i)
973 C if(me.eq.king.or..not.out1file)
974 C & write (iout,*) 'FTORS',ftors
975 C do i=1,ntheta_constr
976 C ii = itheta_constr(i)
977 C thetabound(1,ii) = phi0(i)-drange(i)
978 C thetabound(2,ii) = phi0(i)+drange(i)
980 endif ! ntheta_constr.gt.0
981 endif! with_theta_constr
983 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
984 C write (iout,*) "with_dihed_constr ",with_dihed_constr
989 write (iout,'(a)') 'Boundaries in phi angle sampling:'
991 write (iout,'(a3,i5,2f10.1)')
992 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
998 cd print *,'NNT=',NNT,' NCT=',NCT
999 if (itype(1).eq.ntyp1) nnt=2
1000 if (itype(nres).eq.ntyp1) nct=nct-1
1002 if(me.eq.king.or..not.out1file)
1003 & write (iout,'(a,i3)') 'nsup=',nsup
1005 if (nsup.le.(nct-nnt+1)) then
1006 do i=0,nct-nnt+1-nsup
1007 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1013 & 'Error - sequences to be superposed do not match.'
1016 do i=0,nsup-(nct-nnt+1)
1017 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1019 nstart_sup=nstart_sup+i
1025 & 'Error - sequences to be superposed do not match.'
1028 if (nsup.eq.0) nsup=nct-nnt
1029 if (nstart_sup.eq.0) nstart_sup=nnt
1030 if (nstart_seq.eq.0) nstart_seq=nnt
1031 if(me.eq.king.or..not.out1file)
1032 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1033 & ' nstart_seq=',nstart_seq
1035 c--- Zscore rms -------
1036 if (nz_start.eq.0) nz_start=nnt
1037 if (nz_end.eq.0 .and. nsup.gt.0) then
1039 else if (nz_end.eq.0) then
1042 if(me.eq.king.or..not.out1file)then
1043 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1044 write (iout,*) 'IZ_SC=',iz_sc
1046 c----------------------
1049 if (.not.pdbref) then
1050 call read_angles(inp,*38)
1052 38 write (iout,'(a)') 'Error reading reference structure.'
1054 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1055 stop 'Error reading reference structure'
1057 39 call chainbuild_extconf
1059 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1066 cref(j,i,kkk)=c(j,i)
1069 call contact(.true.,ncont_ref,icont_ref,co)
1073 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1075 if (constr_dist.gt.0) call read_dist_constr
1076 write (iout,*) "After read_dist_constr nhpb",nhpb
1077 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1079 if(me.eq.king.or..not.out1file)
1080 & write (iout,*) 'Contact order:',co
1082 if(me.eq.king.or..not.out1file)
1083 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1086 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1088 if(me.eq.king.or..not.out1file)
1089 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1090 & icont_ref(1,i),' ',
1091 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1095 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1096 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1097 & modecalc.ne.10) then
1098 C If input structure hasn't been supplied from the PDB file read or generate
1100 if (iranconf.eq.0 .and. .not. extconf) then
1101 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1102 & write (iout,'(a)') 'Initial geometry will be read in.'
1104 read(inp,'(8f10.5)',end=36,err=36)
1105 & ((c(l,k),l=1,3),k=1,nres),
1106 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1107 write (iout,*) "Exit READ_CART"
1108 write (iout,'(8f10.5)')
1109 & ((c(l,k),l=1,3),k=1,nres),
1110 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1111 call int_from_cart1(.true.)
1112 write (iout,*) "Finish INT_TO_CART"
1115 dc(j,i)=c(j,i+1)-c(j,i)
1116 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1120 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1122 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1123 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1129 call read_angles(inp,*36)
1130 call chainbuild_extconf
1133 36 write (iout,'(a)') 'Error reading angle file.'
1135 call mpi_finalize( MPI_COMM_WORLD,IERR )
1137 stop 'Error reading angle file.'
1139 else if (extconf) then
1140 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1141 & write (iout,'(a)') 'Extended chain initial geometry.'
1143 theta(i)=90d0*deg2rad
1146 phi(i)=180d0*deg2rad
1149 alph(i)=110d0*deg2rad
1152 omeg(i)=-120d0*deg2rad
1153 if (itype(i).le.0) omeg(i)=-omeg(i)
1155 call chainbuild_extconf
1157 if(me.eq.king.or..not.out1file)
1158 & write (iout,'(a)') 'Random-generated initial geometry.'
1162 if (me.eq.king .or. fg_rank.eq.0 .and. (
1163 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1167 call gen_rand_conf(itmp,*30)
1169 30 write (iout,*) 'Failed to generate random conformation',
1170 & ', itrial=',itrial
1171 write (*,*) 'Processor:',me,
1172 & ' Failed to generate random conformation',
1182 write (iout,'(a,i3,a)') 'Processor:',me,
1183 & ' error in generating random conformation.'
1184 write (*,'(a,i3,a)') 'Processor:',me,
1185 & ' error in generating random conformation.'
1188 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1193 & ' error in generating random conformation.'
1198 elseif (modecalc.eq.4) then
1199 read (inp,'(a)') intinname
1200 open (intin,file=intinname,status='old',err=333)
1201 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1202 & write (iout,'(a)') 'intinname',intinname
1203 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1205 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1207 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1209 stop 'Error opening angle file.'
1213 C Generate distance constraints, if the PDB structure is to be regularized.
1214 if (nthread.gt.0) then
1215 call read_threadbase
1218 if (me.eq.king .or. .not. out1file)
1220 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1221 write (iout,'(/a,i3,a)')
1222 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1223 write (iout,'(20i4)') (iss(i),i=1,ns)
1225 write(iout,*)"Running with dynamic disulfide-bond formation"
1227 write (iout,'(/a/)') 'Pre-formed links are:'
1233 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1234 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1240 if (ns.gt.0.and.dyn_ss) then
1244 forcon(i-nss)=forcon(i)
1251 dyn_ss_mask(iss(i))=.true.
1254 if (i2ndstr.gt.0) call secstrp2dihc
1255 c call geom_to_var(nvar,x)
1256 c call etotal(energia(0))
1257 c call enerprint(energia(0))
1258 c call briefout(0,etot)
1260 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1261 cd write (iout,'(a)') 'Variable list:'
1262 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1264 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1265 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1266 & 'Processor',myrank,': end reading molecular data.'
1271 c--------------------------------------------------------------------------
1272 logical function seq_comp(itypea,itypeb,length)
1274 integer length,itypea(length),itypeb(length)
1277 if (itypea(i).ne.itypeb(i)) then
1285 c-----------------------------------------------------------------------------
1286 subroutine read_bridge
1287 C Read information about disulfide bridges.
1288 implicit real*8 (a-h,o-z)
1289 include 'DIMENSIONS'
1293 include 'COMMON.IOUNITS'
1294 include 'COMMON.GEO'
1295 include 'COMMON.VAR'
1296 include 'COMMON.INTERACT'
1297 include 'COMMON.LOCAL'
1298 include 'COMMON.NAMES'
1299 include 'COMMON.CHAIN'
1300 include 'COMMON.FFIELD'
1301 include 'COMMON.SBRIDGE'
1302 include 'COMMON.HEADER'
1303 include 'COMMON.CONTROL'
1304 include 'COMMON.DBASE'
1305 include 'COMMON.THREAD'
1306 include 'COMMON.TIME1'
1307 include 'COMMON.SETUP'
1308 C Read bridging residues.
1309 read (inp,*) ns,(iss(i),i=1,ns)
1311 if(me.eq.king.or..not.out1file)
1312 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1313 C Check whether the specified bridging residues are cystines.
1315 if (itype(iss(i)).ne.1) then
1316 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1317 & 'Do you REALLY think that the residue ',
1318 & restyp(itype(iss(i))),i,
1319 & ' can form a disulfide bridge?!!!'
1320 write (*,'(2a,i3,a)')
1321 & 'Do you REALLY think that the residue ',
1322 & restyp(itype(iss(i))),i,
1323 & ' can form a disulfide bridge?!!!'
1325 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1330 C Read preformed bridges.
1332 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1334 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1337 C Check if the residues involved in bridges are in the specified list of
1338 C bridging residues.
1341 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1342 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1343 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1344 & ' contains residues present in other pairs.'
1345 write (*,'(a,i3,a)') 'Disulfide pair',i,
1346 & ' contains residues present in other pairs.'
1348 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1354 if (ihpb(i).eq.iss(j)) goto 10
1356 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1359 if (jhpb(i).eq.iss(j)) goto 20
1361 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1367 ihpb(i)=ihpb(i)+nres
1368 jhpb(i)=jhpb(i)+nres
1374 c----------------------------------------------------------------------------
1375 subroutine read_x(kanal,*)
1376 implicit real*8 (a-h,o-z)
1377 include 'DIMENSIONS'
1378 include 'COMMON.GEO'
1379 include 'COMMON.VAR'
1380 include 'COMMON.CHAIN'
1381 include 'COMMON.IOUNITS'
1382 include 'COMMON.CONTROL'
1383 include 'COMMON.LOCAL'
1384 include 'COMMON.INTERACT'
1385 c Read coordinates from input
1387 read(kanal,'(8f10.5)',end=10,err=10)
1388 & ((c(l,k),l=1,3),k=1,nres),
1389 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1392 c(j,2*nres)=c(j,nres)
1394 call int_from_cart1(.false.)
1397 dc(j,i)=c(j,i+1)-c(j,i)
1398 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1402 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1404 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1405 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1413 c----------------------------------------------------------------------------
1414 subroutine read_threadbase
1415 implicit real*8 (a-h,o-z)
1416 include 'DIMENSIONS'
1417 include 'COMMON.IOUNITS'
1418 include 'COMMON.GEO'
1419 include 'COMMON.VAR'
1420 include 'COMMON.INTERACT'
1421 include 'COMMON.LOCAL'
1422 include 'COMMON.NAMES'
1423 include 'COMMON.CHAIN'
1424 include 'COMMON.FFIELD'
1425 include 'COMMON.SBRIDGE'
1426 include 'COMMON.HEADER'
1427 include 'COMMON.CONTROL'
1428 include 'COMMON.DBASE'
1429 include 'COMMON.THREAD'
1430 include 'COMMON.TIME1'
1431 C Read pattern database for threading.
1432 read (icbase,*) nseq
1434 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1435 & nres_base(2,i),nres_base(3,i)
1436 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1438 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1439 c & nres_base(2,i),nres_base(3,i)
1440 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1444 if (weidis.eq.0.0D0) weidis=0.1D0
1453 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1454 write (iout,'(a,i5)') 'nexcl: ',nexcl
1455 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1458 c------------------------------------------------------------------------------
1459 subroutine setup_var
1460 implicit real*8 (a-h,o-z)
1461 include 'DIMENSIONS'
1462 include 'COMMON.IOUNITS'
1463 include 'COMMON.GEO'
1464 include 'COMMON.VAR'
1465 include 'COMMON.INTERACT'
1466 include 'COMMON.LOCAL'
1467 include 'COMMON.NAMES'
1468 include 'COMMON.CHAIN'
1469 include 'COMMON.FFIELD'
1470 include 'COMMON.SBRIDGE'
1471 include 'COMMON.HEADER'
1472 include 'COMMON.CONTROL'
1473 include 'COMMON.DBASE'
1474 include 'COMMON.THREAD'
1475 include 'COMMON.TIME1'
1476 C Set up variable list.
1482 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1484 ialph(i,1)=nvar+nside
1488 if (indphi.gt.0) then
1490 else if (indback.gt.0) then
1495 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1498 c----------------------------------------------------------------------------
1499 subroutine gen_dist_constr
1500 C Generate CA distance constraints.
1501 implicit real*8 (a-h,o-z)
1502 include 'DIMENSIONS'
1503 include 'COMMON.IOUNITS'
1504 include 'COMMON.GEO'
1505 include 'COMMON.VAR'
1506 include 'COMMON.INTERACT'
1507 include 'COMMON.LOCAL'
1508 include 'COMMON.NAMES'
1509 include 'COMMON.CHAIN'
1510 include 'COMMON.FFIELD'
1511 include 'COMMON.SBRIDGE'
1512 include 'COMMON.HEADER'
1513 include 'COMMON.CONTROL'
1514 include 'COMMON.DBASE'
1515 include 'COMMON.THREAD'
1516 include 'COMMON.TIME1'
1517 dimension itype_pdb(maxres)
1518 common /pizda/ itype_pdb
1520 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1521 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1522 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1524 do i=nstart_sup,nstart_sup+nsup-1
1525 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1526 cd & ' seq_pdb', restyp(itype_pdb(i))
1527 do j=i+2,nstart_sup+nsup-1
1529 ihpb(nhpb)=i+nstart_seq-nstart_sup
1530 jhpb(nhpb)=j+nstart_seq-nstart_sup
1532 dhpb(nhpb)=dist(i,j)
1535 cd write (iout,'(a)') 'Distance constraints:'
1540 cd if (ii.gt.nres) then
1545 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1546 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1547 cd & dhpb(i),forcon(i)
1551 c----------------------------------------------------------------------------
1553 implicit real*8 (a-h,o-z)
1554 include 'DIMENSIONS'
1555 include 'COMMON.MAP'
1556 include 'COMMON.IOUNITS'
1557 character*3 angid(4) /'THE','PHI','ALP','OME'/
1558 character*80 mapcard,ucase
1560 read (inp,'(a)') mapcard
1561 mapcard=ucase(mapcard)
1562 if (index(mapcard,'PHI').gt.0) then
1564 else if (index(mapcard,'THE').gt.0) then
1566 else if (index(mapcard,'ALP').gt.0) then
1568 else if (index(mapcard,'OME').gt.0) then
1571 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1572 stop 'Error - illegal variable spec in MAP card.'
1574 call readi (mapcard,'RES1',res1(imap),0)
1575 call readi (mapcard,'RES2',res2(imap),0)
1576 if (res1(imap).eq.0) then
1577 res1(imap)=res2(imap)
1578 else if (res2(imap).eq.0) then
1579 res2(imap)=res1(imap)
1581 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1583 & 'Error - illegal definition of variable group in MAP.'
1584 stop 'Error - illegal definition of variable group in MAP.'
1586 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1587 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1588 call readi(mapcard,'NSTEP',nstep(imap),0)
1589 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1591 & 'Illegal boundary and/or step size specification in MAP.'
1592 stop 'Illegal boundary and/or step size specification in MAP.'
1597 c----------------------------------------------------------------------------
1599 implicit real*8 (a-h,o-z)
1600 include 'DIMENSIONS'
1601 include 'COMMON.IOUNITS'
1602 include 'COMMON.GEO'
1603 include 'COMMON.CSA'
1604 include 'COMMON.BANK'
1605 include 'COMMON.CONTROL'
1607 character*620 mcmcard
1608 call card_concat(mcmcard)
1610 call readi(mcmcard,'NCONF',nconf,50)
1611 call readi(mcmcard,'NADD',nadd,0)
1612 call readi(mcmcard,'JSTART',jstart,1)
1613 call readi(mcmcard,'JEND',jend,1)
1614 call readi(mcmcard,'NSTMAX',nstmax,500000)
1615 call readi(mcmcard,'N0',n0,1)
1616 call readi(mcmcard,'N1',n1,6)
1617 call readi(mcmcard,'N2',n2,4)
1618 call readi(mcmcard,'N3',n3,0)
1619 call readi(mcmcard,'N4',n4,0)
1620 call readi(mcmcard,'N5',n5,0)
1621 call readi(mcmcard,'N6',n6,10)
1622 call readi(mcmcard,'N7',n7,0)
1623 call readi(mcmcard,'N8',n8,0)
1624 call readi(mcmcard,'N9',n9,0)
1625 call readi(mcmcard,'N14',n14,0)
1626 call readi(mcmcard,'N15',n15,0)
1627 call readi(mcmcard,'N16',n16,0)
1628 call readi(mcmcard,'N17',n17,0)
1629 call readi(mcmcard,'N18',n18,0)
1631 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1633 call readi(mcmcard,'NDIFF',ndiff,2)
1634 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1635 call readi(mcmcard,'IS1',is1,1)
1636 call readi(mcmcard,'IS2',is2,8)
1637 call readi(mcmcard,'NRAN0',nran0,4)
1638 call readi(mcmcard,'NRAN1',nran1,2)
1639 call readi(mcmcard,'IRR',irr,1)
1640 call readi(mcmcard,'NSEED',nseed,20)
1641 call readi(mcmcard,'NTOTAL',ntotal,10000)
1642 call reada(mcmcard,'CUT1',cut1,2.0d0)
1643 call reada(mcmcard,'CUT2',cut2,5.0d0)
1644 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1645 call readi(mcmcard,'ICMAX',icmax,3)
1646 call readi(mcmcard,'IRESTART',irestart,0)
1647 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1650 call reada(mcmcard,'DELE',dele,20.0d0)
1651 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1652 call readi(mcmcard,'IREF',iref,0)
1653 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1654 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1655 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1656 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1657 write (iout,*) "NCONF_IN",nconf_in
1660 c----------------------------------------------------------------------------
1661 cfmc subroutine mcmfread
1662 cfmc implicit real*8 (a-h,o-z)
1663 cfmc include 'DIMENSIONS'
1664 cfmc include 'COMMON.MCMF'
1665 cfmc include 'COMMON.IOUNITS'
1666 cfmc include 'COMMON.GEO'
1667 cfmc character*80 ucase
1668 cfmc character*620 mcmcard
1669 cfmc call card_concat(mcmcard)
1671 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1672 cfmc write(iout,*)'MAXRANT=',maxrant
1673 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1674 cfmc write(iout,*)'MAXFAM=',maxfam
1675 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1676 cfmc write(iout,*)'NNET1=',nnet1
1677 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1678 cfmc write(iout,*)'NNET2=',nnet2
1679 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1680 cfmc write(iout,*)'NNET3=',nnet3
1681 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1682 cfmc write(iout,*)'ILASTT=',ilastt
1683 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1684 cfmc write(iout,*)'MAXSTR=',maxstr
1685 cfmc maxstr_f=maxstr/maxfam
1686 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1687 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1688 cfmc write(iout,*)'NMCMF=',nmcmf
1689 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1690 cfmc write(iout,*)'IFOCUS=',ifocus
1691 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1692 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1693 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1694 cfmc write(iout,*)'INTPRT=',intprt
1695 cfmc call readi(mcmcard,'IPRT',iprt,100)
1696 cfmc write(iout,*)'IPRT=',iprt
1697 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1698 cfmc write(iout,*)'IMAXTR=',imaxtr
1699 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1700 cfmc write(iout,*)'MAXEVEN=',maxeven
1701 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1702 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1703 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1704 cfmc write(iout,*)'INIMIN=',inimin
1705 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1706 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1707 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1708 cfmc write(iout,*)'NTHREAD=',nthread
1709 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1710 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1711 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1712 cfmc write(iout,*)'MAXPERT=',maxpert
1713 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1714 cfmc write(iout,*)'IRMSD=',irmsd
1715 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1716 cfmc write(iout,*)'DENEMIN=',denemin
1717 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1718 cfmc write(iout,*)'RCUT1S=',rcut1s
1719 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1720 cfmc write(iout,*)'RCUT1E=',rcut1e
1721 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1722 cfmc write(iout,*)'RCUT2S=',rcut2s
1723 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1724 cfmc write(iout,*)'RCUT2E=',rcut2e
1725 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1726 cfmc write(iout,*)'DPERT1=',d_pert1
1727 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1728 cfmc write(iout,*)'DPERT1A=',d_pert1a
1729 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1730 cfmc write(iout,*)'DPERT2=',d_pert2
1731 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1732 cfmc write(iout,*)'DPERT2A=',d_pert2a
1733 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1734 cfmc write(iout,*)'DPERT2B=',d_pert2b
1735 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1736 cfmc write(iout,*)'DPERT2C=',d_pert2c
1737 cfmc d_pert1=deg2rad*d_pert1
1738 cfmc d_pert1a=deg2rad*d_pert1a
1739 cfmc d_pert2=deg2rad*d_pert2
1740 cfmc d_pert2a=deg2rad*d_pert2a
1741 cfmc d_pert2b=deg2rad*d_pert2b
1742 cfmc d_pert2c=deg2rad*d_pert2c
1743 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1744 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1745 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1746 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1747 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1748 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1749 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1750 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1751 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1752 cfmc write(iout,*)'RCUTINI=',rcutini
1753 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1754 cfmc write(iout,*)'GRAT=',grat
1755 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1756 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1760 c----------------------------------------------------------------------------
1762 implicit real*8 (a-h,o-z)
1763 include 'DIMENSIONS'
1764 include 'COMMON.MCM'
1765 include 'COMMON.MCE'
1766 include 'COMMON.IOUNITS'
1768 character*320 mcmcard
1769 call card_concat(mcmcard)
1770 call readi(mcmcard,'MAXACC',maxacc,100)
1771 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1772 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1773 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1774 call readi(mcmcard,'MAXREPM',maxrepm,200)
1775 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1776 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1777 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1778 call reada(mcmcard,'E_UP',e_up,5.0D0)
1779 call reada(mcmcard,'DELTE',delte,0.1D0)
1780 call readi(mcmcard,'NSWEEP',nsweep,5)
1781 call readi(mcmcard,'NSTEPH',nsteph,0)
1782 call readi(mcmcard,'NSTEPC',nstepc,0)
1783 call reada(mcmcard,'TMIN',tmin,298.0D0)
1784 call reada(mcmcard,'TMAX',tmax,298.0D0)
1785 call readi(mcmcard,'NWINDOW',nwindow,0)
1786 call readi(mcmcard,'PRINT_MC',print_mc,0)
1787 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1788 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1789 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1790 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1791 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1792 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1793 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1794 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1795 if (nwindow.gt.0) then
1796 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1798 winlen(i)=winend(i)-winstart(i)+1
1801 if (tmax.lt.tmin) tmax=tmin
1802 if (tmax.eq.tmin) then
1806 if (nstepc.gt.0 .and. nsteph.gt.0) then
1807 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1808 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1810 C Probabilities of different move types
1811 sumpro_type(0)=0.0D0
1812 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1813 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1814 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1815 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1816 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1817 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1818 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1820 print *,'i',i,' sumprotype',sumpro_type(i)
1821 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1822 print *,'i',i,' sumprotype',sumpro_type(i)
1826 c----------------------------------------------------------------------------
1827 subroutine read_minim
1828 implicit real*8 (a-h,o-z)
1829 include 'DIMENSIONS'
1830 include 'COMMON.MINIM'
1831 include 'COMMON.IOUNITS'
1833 character*320 minimcard
1834 call card_concat(minimcard)
1835 call readi(minimcard,'MAXMIN',maxmin,2000)
1836 call readi(minimcard,'MAXFUN',maxfun,5000)
1837 call readi(minimcard,'MINMIN',minmin,maxmin)
1838 call readi(minimcard,'MINFUN',minfun,maxmin)
1839 call reada(minimcard,'TOLF',tolf,1.0D-2)
1840 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1841 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1842 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1843 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1844 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1845 & 'Options in energy minimization:'
1846 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1847 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1848 & 'MinMin:',MinMin,' MinFun:',MinFun,
1849 & ' TolF:',TolF,' RTolF:',RTolF
1852 c----------------------------------------------------------------------------
1853 subroutine read_angles(kanal,*)
1854 implicit real*8 (a-h,o-z)
1855 include 'DIMENSIONS'
1856 include 'COMMON.GEO'
1857 include 'COMMON.VAR'
1858 include 'COMMON.CHAIN'
1859 include 'COMMON.IOUNITS'
1860 include 'COMMON.CONTROL'
1861 c Read angles from input
1863 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1864 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1865 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1866 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1869 c 9/7/01 avoid 180 deg valence angle
1870 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1872 theta(i)=deg2rad*theta(i)
1873 phi(i)=deg2rad*phi(i)
1874 alph(i)=deg2rad*alph(i)
1875 omeg(i)=deg2rad*omeg(i)
1880 c----------------------------------------------------------------------------
1881 subroutine reada(rekord,lancuch,wartosc,default)
1883 character*(*) rekord,lancuch
1884 double precision wartosc,default
1887 iread=index(rekord,lancuch)
1888 if (iread.eq.0) then
1892 iread=iread+ilen(lancuch)+1
1893 read (rekord(iread:),*,err=10,end=10) wartosc
1898 c----------------------------------------------------------------------------
1899 subroutine readi(rekord,lancuch,wartosc,default)
1901 character*(*) rekord,lancuch
1902 integer wartosc,default
1905 iread=index(rekord,lancuch)
1906 if (iread.eq.0) then
1910 iread=iread+ilen(lancuch)+1
1911 read (rekord(iread:),*,err=10,end=10) wartosc
1916 c----------------------------------------------------------------------------
1917 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1920 integer tablica(dim),default
1921 character*(*) rekord,lancuch
1928 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1929 if (iread.eq.0) return
1930 iread=iread+ilen(lancuch)+1
1931 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1934 c----------------------------------------------------------------------------
1935 subroutine multreada(rekord,lancuch,tablica,dim,default)
1938 double precision tablica(dim),default
1939 character*(*) rekord,lancuch
1946 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1947 if (iread.eq.0) return
1948 iread=iread+ilen(lancuch)+1
1949 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1952 c----------------------------------------------------------------------------
1953 subroutine openunits
1954 implicit real*8 (a-h,o-z)
1955 include 'DIMENSIONS'
1958 character*16 form,nodename
1961 include 'COMMON.SETUP'
1962 include 'COMMON.IOUNITS'
1964 include 'COMMON.CONTROL'
1965 integer lenpre,lenpot,ilen,lentmp
1967 character*3 out1file_text,ucase
1970 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1971 call getenv_loc("PREFIX",prefix)
1973 call getenv_loc("POT",pot)
1974 call getenv_loc("DIRTMP",tmpdir)
1975 call getenv_loc("CURDIR",curdir)
1976 call getenv_loc("OUT1FILE",out1file_text)
1977 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1978 out1file_text=ucase(out1file_text)
1979 if (out1file_text(1:1).eq."Y") then
1982 out1file=fg_rank.gt.0
1987 if (lentmp.gt.0) then
1988 write (*,'(80(1h!))')
1989 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1990 write (*,'(80(1h!))')
1991 write (*,*)"All output files will be on node /tmp directory."
1993 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1994 if (me.eq.king) then
1995 write (*,*) "The master node is ",nodename
1996 else if (fg_rank.eq.0) then
1997 write (*,*) "I am the CG slave node ",nodename
1999 write (*,*) "I am the FG slave node ",nodename
2002 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2003 lenpre = lentmp+lenpre+1
2005 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2006 C Get the names and open the input files
2007 #if defined(WINIFL) || defined(WINPGI)
2008 open(1,file=pref_orig(:ilen(pref_orig))//
2009 & '.inp',status='old',readonly,shared)
2010 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2011 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2012 C Get parameter filenames and open the parameter files.
2013 call getenv_loc('BONDPAR',bondname)
2014 open (ibond,file=bondname,status='old',readonly,shared)
2015 call getenv_loc('THETPAR',thetname)
2016 open (ithep,file=thetname,status='old',readonly,shared)
2017 call getenv_loc('ROTPAR',rotname)
2018 open (irotam,file=rotname,status='old',readonly,shared)
2019 call getenv_loc('TORPAR',torname)
2020 open (itorp,file=torname,status='old',readonly,shared)
2021 call getenv_loc('TORDPAR',tordname)
2022 open (itordp,file=tordname,status='old',readonly,shared)
2023 call getenv_loc('TORKCC',torkccname)
2024 open (itorkcc,file=torkccname,status='old',readonly,shared)
2025 call getenv_loc('THETKCC',thetkccname)
2026 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2027 call getenv_loc('FOURIER',fouriername)
2028 open (ifourier,file=fouriername,status='old',readonly,shared)
2029 call getenv_loc('ELEPAR',elename)
2030 open (ielep,file=elename,status='old',readonly,shared)
2031 call getenv_loc('SIDEPAR',sidename)
2032 open (isidep,file=sidename,status='old',readonly,shared)
2033 call getenv_loc('LIPTRANPAR',liptranname)
2034 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2035 #elif (defined CRAY) || (defined AIX)
2036 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2038 c print *,"Processor",myrank," opened file 1"
2039 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2040 c print *,"Processor",myrank," opened file 9"
2041 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2042 C Get parameter filenames and open the parameter files.
2043 call getenv_loc('BONDPAR',bondname)
2044 open (ibond,file=bondname,status='old',action='read')
2045 c print *,"Processor",myrank," opened file IBOND"
2046 call getenv_loc('THETPAR',thetname)
2047 open (ithep,file=thetname,status='old',action='read')
2048 c print *,"Processor",myrank," opened file ITHEP"
2049 call getenv_loc('ROTPAR',rotname)
2050 open (irotam,file=rotname,status='old',action='read')
2051 c print *,"Processor",myrank," opened file IROTAM"
2052 call getenv_loc('TORPAR',torname)
2053 open (itorp,file=torname,status='old',action='read')
2054 c print *,"Processor",myrank," opened file ITORP"
2055 call getenv_loc('TORDPAR',tordname)
2056 open (itordp,file=tordname,status='old',action='read')
2057 call getenv_loc('TORKCC',torkccname)
2058 open (itorkcc,file=torkccname,status='old',action='read')
2059 call getenv_loc('THETKCC',thetkccname)
2060 open (ithetkcc,file=thetkccname,status='old',action='read')
2061 c print *,"Processor",myrank," opened file ITORDP"
2062 call getenv_loc('SCCORPAR',sccorname)
2063 open (isccor,file=sccorname,status='old',action='read')
2064 c print *,"Processor",myrank," opened file ISCCOR"
2065 call getenv_loc('FOURIER',fouriername)
2066 open (ifourier,file=fouriername,status='old',action='read')
2067 c print *,"Processor",myrank," opened file IFOURIER"
2068 call getenv_loc('ELEPAR',elename)
2069 open (ielep,file=elename,status='old',action='read')
2070 c print *,"Processor",myrank," opened file IELEP"
2071 call getenv_loc('SIDEPAR',sidename)
2072 open (isidep,file=sidename,status='old',action='read')
2073 call getenv_loc('LIPTRANPAR',liptranname)
2074 open (iliptranpar,file=liptranname,status='old',action='read')
2075 c print *,"Processor",myrank," opened file ISIDEP"
2076 c print *,"Processor",myrank," opened parameter files"
2078 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2079 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2080 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2081 C Get parameter filenames and open the parameter files.
2082 call getenv_loc('BONDPAR',bondname)
2083 open (ibond,file=bondname,status='old')
2084 call getenv_loc('THETPAR',thetname)
2085 open (ithep,file=thetname,status='old')
2086 call getenv_loc('ROTPAR',rotname)
2087 open (irotam,file=rotname,status='old')
2088 call getenv_loc('TORPAR',torname)
2089 open (itorp,file=torname,status='old')
2090 call getenv_loc('TORDPAR',tordname)
2091 open (itordp,file=tordname,status='old')
2092 call getenv_loc('TORKCC',torkccname)
2093 open (itorkcc,file=torkccname,status='old')
2094 call getenv_loc('THETKCC',thetkccname)
2095 open (ithetkcc,file=thetkccname,status='old')
2096 call getenv_loc('SCCORPAR',sccorname)
2097 open (isccor,file=sccorname,status='old')
2098 call getenv_loc('FOURIER',fouriername)
2099 open (ifourier,file=fouriername,status='old')
2100 call getenv_loc('ELEPAR',elename)
2101 open (ielep,file=elename,status='old')
2102 call getenv_loc('SIDEPAR',sidename)
2103 open (isidep,file=sidename,status='old')
2104 call getenv_loc('LIPTRANPAR',liptranname)
2105 open (iliptranpar,file=liptranname,status='old')
2107 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2109 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2110 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2111 C Get parameter filenames and open the parameter files.
2112 call getenv_loc('BONDPAR',bondname)
2113 open (ibond,file=bondname,status='old',readonly)
2114 call getenv_loc('THETPAR',thetname)
2115 open (ithep,file=thetname,status='old',readonly)
2116 call getenv_loc('ROTPAR',rotname)
2117 open (irotam,file=rotname,status='old',readonly)
2118 call getenv_loc('TORPAR',torname)
2119 open (itorp,file=torname,status='old',readonly)
2120 call getenv_loc('TORDPAR',tordname)
2121 open (itordp,file=tordname,status='old',readonly)
2122 call getenv_loc('TORKCC',torkccname)
2123 open (itorkcc,file=torkccname,status='old',readonly)
2124 call getenv_loc('THETKCC',thetkccname)
2125 open (ithetkcc,file=thetkccname,status='old',readonly)
2126 call getenv_loc('SCCORPAR',sccorname)
2127 open (isccor,file=sccorname,status='old',readonly)
2129 call getenv_loc('THETPARPDB',thetname_pdb)
2130 print *,"thetname_pdb ",thetname_pdb
2131 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2132 print *,ithep_pdb," opened"
2134 call getenv_loc('FOURIER',fouriername)
2135 open (ifourier,file=fouriername,status='old',readonly)
2136 call getenv_loc('ELEPAR',elename)
2137 open (ielep,file=elename,status='old',readonly)
2138 call getenv_loc('SIDEPAR',sidename)
2139 open (isidep,file=sidename,status='old',readonly)
2140 call getenv_loc('LIPTRANPAR',liptranname)
2141 open (iliptranpar,file=liptranname,status='old',action='read')
2143 call getenv_loc('ROTPARPDB',rotname_pdb)
2144 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2147 call getenv_loc('TUBEPAR',tubename)
2148 open (itube,file=tubename,status='old',readonly)
2151 C 8/9/01 In the newest version SCp interaction constants are read from a file
2152 C Use -DOLDSCP to use hard-coded constants instead.
2154 call getenv_loc('SCPPAR',scpname)
2155 #if defined(WINIFL) || defined(WINPGI)
2156 open (iscpp,file=scpname,status='old',readonly,shared)
2157 #elif (defined CRAY) || (defined AIX)
2158 open (iscpp,file=scpname,status='old',action='read')
2160 open (iscpp,file=scpname,status='old')
2162 open (iscpp,file=scpname,status='old',readonly)
2165 call getenv_loc('PATTERN',patname)
2166 #if defined(WINIFL) || defined(WINPGI)
2167 open (icbase,file=patname,status='old',readonly,shared)
2168 #elif (defined CRAY) || (defined AIX)
2169 open (icbase,file=patname,status='old',action='read')
2171 open (icbase,file=patname,status='old')
2173 open (icbase,file=patname,status='old',readonly)
2176 C Open output file only for CG processes
2177 c print *,"Processor",myrank," fg_rank",fg_rank
2178 if (fg_rank.eq.0) then
2180 if (nodes.eq.1) then
2183 npos = dlog10(dfloat(nodes-1))+1
2185 if (npos.lt.3) npos=3
2186 write (liczba,'(i1)') npos
2187 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2189 write (liczba,form) me
2190 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2191 & liczba(:ilen(liczba))
2192 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2194 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2196 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2197 & liczba(:ilen(liczba))//'.mol2'
2198 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2199 & liczba(:ilen(liczba))//'.stat'
2201 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2202 & //liczba(:ilen(liczba))//'.stat')
2203 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2206 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2207 & liczba(:ilen(liczba))//'.const'
2212 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2213 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2214 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2215 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2216 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2218 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2220 rest2name=prefix(:ilen(prefix))//'.rst'
2222 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2225 #if defined(AIX) || defined(PGI)
2226 if (me.eq.king .or. .not. out1file)
2227 & open(iout,file=outname,status='unknown')
2229 if (fg_rank.gt.0) then
2230 write (liczba,'(i3.3)') myrank/nfgtasks
2231 write (ll,'(bz,i3.3)') fg_rank
2232 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2237 open(igeom,file=intname,status='unknown',position='append')
2238 open(ipdb,file=pdbname,status='unknown')
2239 open(imol2,file=mol2name,status='unknown')
2240 open(istat,file=statname,status='unknown',position='append')
2242 c1out open(iout,file=outname,status='unknown')
2245 if (me.eq.king .or. .not.out1file)
2246 & open(iout,file=outname,status='unknown')
2248 if (fg_rank.gt.0) then
2249 write (liczba,'(i3.3)') myrank/nfgtasks
2250 write (ll,'(bz,i3.3)') fg_rank
2251 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2256 open(igeom,file=intname,status='unknown',access='append')
2257 open(ipdb,file=pdbname,status='unknown')
2258 open(imol2,file=mol2name,status='unknown')
2259 open(istat,file=statname,status='unknown',access='append')
2261 c1out open(iout,file=outname,status='unknown')
2264 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2265 csa_seed=prefix(:lenpre)//'.CSA.seed'
2266 csa_history=prefix(:lenpre)//'.CSA.history'
2267 csa_bank=prefix(:lenpre)//'.CSA.bank'
2268 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2269 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2270 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2271 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2272 csa_int=prefix(:lenpre)//'.int'
2273 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2274 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2275 csa_in=prefix(:lenpre)//'.CSA.in'
2276 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2279 write (iout,'(80(1h-))')
2280 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2281 write (iout,'(80(1h-))')
2282 write (iout,*) "Input file : ",
2283 & pref_orig(:ilen(pref_orig))//'.inp'
2284 write (iout,*) "Output file : ",
2285 & outname(:ilen(outname))
2287 write (iout,*) "Sidechain potential file : ",
2288 & sidename(:ilen(sidename))
2290 write (iout,*) "SCp potential file : ",
2291 & scpname(:ilen(scpname))
2293 write (iout,*) "Electrostatic potential file : ",
2294 & elename(:ilen(elename))
2295 write (iout,*) "Cumulant coefficient file : ",
2296 & fouriername(:ilen(fouriername))
2297 write (iout,*) "Torsional parameter file : ",
2298 & torname(:ilen(torname))
2299 write (iout,*) "Double torsional parameter file : ",
2300 & tordname(:ilen(tordname))
2301 write (iout,*) "SCCOR parameter file : ",
2302 & sccorname(:ilen(sccorname))
2303 write (iout,*) "Bond & inertia constant file : ",
2304 & bondname(:ilen(bondname))
2305 write (iout,*) "Bending parameter file : ",
2306 & thetname(:ilen(thetname))
2307 write (iout,*) "Rotamer parameter file : ",
2308 & rotname(:ilen(rotname))
2309 write (iout,*) "Thetpdb parameter file : ",
2310 & thetname_pdb(:ilen(thetname_pdb))
2311 write (iout,*) "Threading database : ",
2312 & patname(:ilen(patname))
2314 &write (iout,*)" DIRTMP : ",
2316 write (iout,'(80(1h-))')
2320 c----------------------------------------------------------------------------
2321 subroutine card_concat(card)
2322 implicit real*8 (a-h,o-z)
2323 include 'DIMENSIONS'
2324 include 'COMMON.IOUNITS'
2326 character*80 karta,ucase
2328 read (inp,'(a)') karta
2331 do while (karta(80:80).eq.'&')
2332 card=card(:ilen(card)+1)//karta(:79)
2333 read (inp,'(a)') karta
2336 card=card(:ilen(card)+1)//karta
2339 c----------------------------------------------------------------------------------
2341 implicit real*8 (a-h,o-z)
2342 include 'DIMENSIONS'
2343 include 'COMMON.CHAIN'
2344 include 'COMMON.IOUNITS'
2346 open(irest2,file=rest2name,status='unknown')
2347 read(irest2,*) totT,EK,potE,totE,t_bath
2350 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2353 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2356 read (irest2,*) iset
2361 c---------------------------------------------------------------------------------
2362 subroutine read_fragments
2363 implicit real*8 (a-h,o-z)
2364 include 'DIMENSIONS'
2368 include 'COMMON.SETUP'
2369 include 'COMMON.CHAIN'
2370 include 'COMMON.IOUNITS'
2372 include 'COMMON.CONTROL'
2373 read(inp,*) nset,nfrag,npair,nfrag_back
2374 if(me.eq.king.or..not.out1file)
2375 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2376 & " nfrag_back",nfrag_back
2378 read(inp,*) mset(iset)
2380 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2382 if(me.eq.king.or..not.out1file)
2383 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2384 & ifrag(2,i,iset), qinfrag(i,iset)
2387 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2389 if(me.eq.king.or..not.out1file)
2390 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2391 & ipair(2,i,iset), qinpair(i,iset)
2394 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2395 & wfrag_back(3,i,iset),
2396 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2397 if(me.eq.king.or..not.out1file)
2398 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2399 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2404 C---------------------------------------------------------------------------
2405 subroutine read_afminp
2406 implicit real*8 (a-h,o-z)
2407 include 'DIMENSIONS'
2411 include 'COMMON.SETUP'
2412 include 'COMMON.CONTROL'
2413 include 'COMMON.CHAIN'
2414 include 'COMMON.IOUNITS'
2415 include 'COMMON.SBRIDGE'
2416 character*320 afmcard
2418 call card_concat(afmcard)
2419 call readi(afmcard,"BEG",afmbeg,0)
2420 call readi(afmcard,"END",afmend,0)
2421 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2422 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2423 print *,'FORCE=' ,forceAFMconst
2424 CCCC NOW PROPERTIES FOR AFM
2427 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2429 distafminit=dsqrt(distafminit)
2430 print *,'initdist',distafminit
2434 c-------------------------------------------------------------------------------
2435 subroutine read_dist_constr
2436 implicit real*8 (a-h,o-z)
2437 include 'DIMENSIONS'
2441 include 'COMMON.SETUP'
2442 include 'COMMON.CONTROL'
2443 include 'COMMON.CHAIN'
2444 include 'COMMON.IOUNITS'
2445 include 'COMMON.SBRIDGE'
2446 integer ifrag_(2,100),ipair_(2,100)
2447 double precision wfrag_(100),wpair_(100)
2448 character*500 controlcard
2450 write (iout,*) "Calling read_dist_constr"
2451 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2453 call card_concat(controlcard)
2454 call readi(controlcard,"NFRAG",nfrag_,0)
2455 call readi(controlcard,"NPAIR",npair_,0)
2456 call readi(controlcard,"NDIST",ndist_,0)
2457 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2458 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2459 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2460 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2461 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2462 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2463 c write (iout,*) "IFRAG"
2465 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2467 c write (iout,*) "IPAIR"
2469 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2473 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2474 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2475 & ifrag_(2,i)=nstart_sup+nsup-1
2476 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2478 if (wfrag_(i).gt.0.0d0) then
2479 do j=ifrag_(1,i),ifrag_(2,i)-1
2480 do k=j+1,ifrag_(2,i)
2481 c write (iout,*) "j",j," k",k
2483 if (constr_dist.eq.1) then
2488 forcon(nhpb)=wfrag_(i)
2489 else if (constr_dist.eq.2) then
2490 if (ddjk.le.dist_cut) then
2495 forcon(nhpb)=wfrag_(i)
2502 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2505 if (.not.out1file .or. me.eq.king)
2506 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2507 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2509 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2510 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2517 if (wpair_(i).gt.0.0d0) then
2525 do j=ifrag_(1,ii),ifrag_(2,ii)
2526 do k=ifrag_(1,jj),ifrag_(2,jj)
2530 forcon(nhpb)=wpair_(i)
2531 dhpb(nhpb)=dist(j,k)
2533 if (.not.out1file .or. me.eq.king)
2534 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2535 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2537 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2538 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2546 if (constr_dist.eq.11) then
2547 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2548 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2549 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2552 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2553 & ibecarb(i),forcon(nhpb+1)
2555 if (forcon(nhpb+1).gt.0.0d0) then
2557 if (ibecarb(i).gt.0) then
2558 ihpb(i)=ihpb(i)+nres
2559 jhpb(i)=jhpb(i)+nres
2561 if (dhpb(nhpb).eq.0.0d0)
2562 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2564 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2565 C if (forcon(nhpb+1).gt.0.0d0) then
2567 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2569 if (.not.out1file .or. me.eq.king)
2570 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2571 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2573 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2574 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2581 c-------------------------------------------------------------------------------
2583 subroutine flush(iu)
2588 subroutine flush(iu)
2593 c------------------------------------------------------------------------------
2594 subroutine copy_to_tmp(source)
2595 include "DIMENSIONS"
2596 include "COMMON.IOUNITS"
2597 character*(*) source
2598 character* 256 tmpfile
2602 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2603 inquire(file=tmpfile,exist=ex)
2605 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2606 & " to temporary directory..."
2607 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2608 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2612 c------------------------------------------------------------------------------
2613 subroutine move_from_tmp(source)
2614 include "DIMENSIONS"
2615 include "COMMON.IOUNITS"
2616 character*(*) source
2619 write (*,*) "Moving ",source(:ilen(source)),
2620 & " from temporary directory to working directory"
2621 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2622 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2625 c------------------------------------------------------------------------------
2626 subroutine random_init(seed)
2628 C Initialize random number generator
2630 implicit real*8 (a-h,o-z)
2631 include 'DIMENSIONS'
2634 logical OKRandom, prng_restart
2636 integer iseed_array(4)
2638 include 'COMMON.IOUNITS'
2639 include 'COMMON.TIME1'
2640 include 'COMMON.THREAD'
2641 include 'COMMON.SBRIDGE'
2642 include 'COMMON.CONTROL'
2643 include 'COMMON.MCM'
2644 include 'COMMON.MAP'
2645 include 'COMMON.HEADER'
2646 include 'COMMON.CSA'
2647 include 'COMMON.CHAIN'
2648 include 'COMMON.MUCA'
2650 include 'COMMON.FFIELD'
2651 include 'COMMON.SETUP'
2652 iseed=-dint(dabs(seed))
2653 if (iseed.eq.0) then
2654 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2655 & 'Random seed undefined. The program will stop.'
2656 write (*,'(/80(1h*)/20x,a/80(1h*))')
2657 & 'Random seed undefined. The program will stop.'
2659 call mpi_finalize(mpi_comm_world,ierr)
2661 stop 'Bad random seed.'
2664 if (fg_rank.eq.0) then
2668 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2669 OKRandom = prng_restart(me,iseed)
2672 tmp=65536.0d0**(4-i)
2673 iseed_array(i) = dint(seed/tmp)
2674 seed=seed-iseed_array(i)*tmp
2677 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2678 & (iseed_array(i),i=1,4)
2679 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2680 & (iseed_array(i),i=1,4)
2681 OKRandom = prng_restart(me,iseed_array)
2684 c r1 = prng_next(me)
2685 r1=ran_number(0.0D0,1.0D0)
2687 & write (iout,*) 'ran_num',r1
2688 if (r1.lt.0.0d0) OKRandom=.false.
2690 if (.not.OKRandom) then
2691 write (iout,*) 'PRNG IS NOT WORKING!!!'
2692 print *,'PRNG IS NOT WORKING!!!'
2695 call mpi_abort(mpi_comm_world,error_msg,ierr)
2698 write (iout,*) 'too many processors for parallel prng'
2699 write (*,*) 'too many processors for parallel prng'
2707 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)